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1. Introduction. It has been shown by J. Schwinger that a special class of bound- 
ary value problems in electrodynamics can be formulated mathematically as Wiener- 
Hopf‘ integral equations. These problems may be described as follows. A plane wave 
is incident upon a number of semi-infinite parallel metallic structures of zero thickness 
and perfect conductivity. By parallel structures we mean parallel planes or cylinders 
with parallel axes. It is then possible to express the electric or magnetic field at all 
points in space in terms of the surface current density on the metal with the aid of an 
appropriate Green’s function. The vanishing of the components of the electric field 
which are tangential to the semi-infinite cylindrical metallic surfaces, leads to a sys- 
tem of inhomogeneous integral equations for the various surface current densities. 
This system of integral equations assumes the general form 


g(x) = D0 Kix —yfdydy, «>0, i=1,---,K, 
j=-1/4 0 
where the f;(y) are unknown functions, while the K;;(x) and g,;(x) are known. The 
particular problem which we shall discuss below possesses certain periodicities, and 
for this case we find it possible to reduce the system to a single integral equation of 
the form 


ow 


s(z) = f K(x— yfy)dy, 2 >0, (1.1) 


that is, an inhomogeneous Wiener-Hopf integral equation. Here f(y) is unknown, 


while K(x) and g(x) are known functions. 
The advantage of formulating this particular class of boundary value problems as 
Wiener-Hopf integral equations is that such equations are susceptible to a rigorous 
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solution. We may thus find the functional form of the various surface current densi- 
ties as well as the electric field. However, in such problems as we have described above, 
the physically interesting quantities may be calculated from the far field and these 
quantities in turn are closely related to the Fourier transform of the surface current 
densities. Since Eq. (1.1) is solved by transform techniques, these quantities can be 
obtained immediately. 

The problem which we treat here is the following. A plane monochromatic electro- 
magnetic wave whose direction of propagation lies in the plane of the paper, is inci- 
dent upon an infinite set of staggered, equally spaced, semi-infinite metallic plates 
of zero thickness and perfect conductivity. These plates extend indefinitely in a direc- 
tion perpendicular to the plane of the paper. (See Fig. 1 for a side view.) The angle of 
stagger with respect to a fixed direction (that of the cross section of the plates in 
Fig. 1) is a, while the direction of propagation with respect to this fixed line is 0, 
where a—7 <6 <aand 0<aSz/2. This structure has some properties which are analo- 
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gous to those of metal mirrors and gratings. Thus when it is excited by a plane wave 
with arbitrary direction of propagation, there will be reflected plane waves in certain 
directions depending on the relative dimensions, the wave length and the direction 
of incidence. 

2. Formulation of the problem. We assume that the electric field of the incident 
wave has only one component, namely, the component which is perpendicular to the 
plane of the paper. Since the incident electric field is independent of y and the bound- 
ary conditions on the plates must be fulfilled independently of y, no other components 
of the electric field will be excited. Thus all components of the magnetic field can 
be derived from this single component of the electric field E,(x, 2) =@(x, 2). For this 
case both of the components of the magnetic field lie in the plane of the paper and 
we shall refer to this problem as an “#7 plane” problem. 

If we now write the Maxwell equations’ in the form 

Vv X E = ikH 
and 
VXH = —-itkE, 
where k = 27/X, and ) is the free space wave-length, we see immediately that the only 
components of the magnetic field are 

5 The time dependence of all field quantities is taken to be e~**«t and may therefore be suppressed. 
c is the velocity of light. In the engineering literature, the time dependence is written as exp(ikct). In 
order to convert our final results to standard engineering form, one merely replaces i by —j. 
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0 

tkH, = — = 

and y 
0 

ikH, = ~ 

Ox 


Upon eliminating 7, and H, from the above equations we obtain the two dimensional 
wave equation, 

ee 

Ox? 02? 
which is to be solved subject to the boundary condition, ¢=0 on the metal plates 
since @ is the tangential component of the electric field. There are also conditions at 
infinity on the function ¢(x, z) which we shall discuss later when we have need of 
them. 

We now formulate the equation which expresses the electric field in terms of the 
surface current density on the metal plates. To this end, we start by modifying the 
structure in Fig. 1, so that there are now only a finite number of parallel plates, each 
of which is taken to be finite in length. The length of each plate is such that the ampli- 
tudes of the attenuated modes are negligibly small relative to the amplitude of the 
propagated mode in the parallel plate region before the end of the structure is 
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reached. (See Fig. 2 for a side view.) If we employ the free space Green’s function, we 
may express $(x, z) in terms of 0¢/dn, the normal derivative on the metallic plates. 


We have from Green’s theorem 


Lots 
x, 2) = — ¢— js’, 
O(a c\ On’ oon’ 
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where the contour C is the one indicated by the dotted line in Fig. 2, ds’ is the element 
of arc length along it and G(x, z, x’, 2’) is the free space Green’s function. The outer 
boundary of the contour C is taken to be a circle of large radius. This is merely for 
convenience and the outer boundary might have been any other closed curve. 
G(x, 2, x’, 2’) satisfies the homogeneous wave equation 


0G dG 





save for the point x =x’, z=2’. At this point 


© 9G |z=2'+0 
f — | dz’ = —1 
—o Ox awe 
and 
oe 0G |2=2’+0 
f —| dx’ = — 1. 
— Oz lge2’—0 


This may be expressed symbolically by saying that G(x, 2, x’, 2’) satisfies the inhomo- 
geneous wave equation 

VG aG 
rs 


+ RG = — i(x — x’)i(z — 2’), 
Ox? 02” ( m= 4) 








where 5(x—x’) is the Dirac delta function and is zero everywhere save at x =x’, where 
it becomes infinite in such a fashion as to make the integral 


f 6(x — x’)dx’ = 1. 


On the plates ¢(x, 2) =0, while 0¢/0n’ is the tangential component of the magnetic 
field on the plates. Since the tangential component of the magnetic field suffers a dis- 
continuity which is proportional to the surface current density when we go from one 
side of a given plate to the other side of it, the only contribution we get from the in- 


tegration along the metallic plates is 


qa 
z. f G(x, z, ma, 2')I m(z’)d2’, 
m=p 


and the limits of integration are those which cover the full length of each plate. The 
sum is carried out over the finite number of plates as shown in Fig. 2. I,,(z) is propor- 
tional to the surface current density on the mth metal plate. There is complete can- 
cellation of the integrals taken along the paths which lead from one plate to the next 
or which lead from the end plates to the large circle enclosing all of the plates. 

We now calculate the contribution from the large circle. In the first place, the free 
space Green’s function which represents an outgoing wave for /x?+-22>>V/x"+2” 
is G(x, 2, x’, 2!) = (4/4) HY [kV/(x—x’)?+ (2—2’)?] where H® is the Hankel function 
of the first kind. The contribution from the large circle is 


d¢(r’, B’) aG(r, r’, B, B’) 
f | Ger, rt. #)— — $(7’, B’) ee ae (2.1) 
0 Or’ Or 
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where 





(1) 


G(r, 7’, B, B’) = = HS [k\/7® + 7? — 297’ cos (B— BY] 


and x =r sin 8B, =r cos B. If we now expand G(r, r’, 8, 8’) in terms of cylindrical waves 
we have 
eo 


G(r, 1’, B, B’) = Lo AM (ker )In(krjeim PO), or <r’, 


mMm=—o0 


Furthermore, for any point outside of the region of the plates 
(x, 2) = cik(s cos O+zsin 6) 4 im anH( kr)ein®, (2.2) 


where the first term represents the incident plane wave whose direction of propaga- 
tion is 6, while the second term represents the scattered wave. We shall not be inter- 
ested in the explicit form of the a,’s and indeed, we shall show that they do not enter 
explicitly into the formulation of the integral equation. The expression for the plane 
wave, e** ls cos +z sin 6) may be expanded in terms of cylindrical waves by noting that 


co) 
eik(2 cos O+2 sin 0) — eikr cos (0-B)" = % ) retimnl2 7. (kr)eim(o-8), 
m=—o0 
If we now evaluate the integrals in (2.1) we get immediately 
eik(z cos 6+zsin 6) — Pinel x, 2) 


i.e., the incident field. 
For our final equation we then have 


(x, 2) = dine(x, 2) + > f I n(2’)G(x, 2, ma, 2’)dz’. 


If we now let g become positively infinite, p negatively infinite, and let each plate 
extend indefinitely to the right, we can then express $(x, 2), the y component of the 
electric field, in terms of the incident field and the surface current density on’ the 


plates, that is, 


$(x, 2) = Pine(X, 2) + oe >>3 In(2’)H [/(3 — 2’)? + (x — ma)?|dz’, (2.3) 


m=-—-oO m b 4 





where a=) tan a. We now impose the electromagnetic boundary condition, namely 
that (x, 2) vanishes on the metallic plates, and we get a system of simultaneous in- 
tegral equations of the Wiener-Hopf type for J,,(z). That is, for x =na 





1 co) [o} 
0 = ¢dine(na, z) +— >> In(2’)H® [k/(z — 2)? + (n — m)?a*]dz’ (2.4) 
m=—o mb 
for all m with s>nb, n=0, +1, +2, - - - 6 Due to the periodic nature of the structure, 
the infinite set of simultaneous integral equations can be cast into the form (1.1). 





® It is possible to obtain the integral equation (2.3) directly from the infinite structure indicated in 
Fig. 1. We have intentionally avoided this because it requires a more detailed knowledge of the field at 


infinity. 
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We close our discussion of the formulation of the integral equation (2.3) with some 
remarks about the range of values of a/X which is allowed. In the parallel plate regions 
for z large and positive, ¢(x, z) is asymptotic to sin (1x/a)e** where x =+/k?— (/a)?. 
If now k<7/a, i.e., \>2a, x will be pure imaginary and hence for z sufficiently 
large and positive, $(x, 2) will vanish exponentially. In this case, the parallel plate 
regions cannot sustain a propagating mode. If k>7/a, i.e., 2a <X, then x is real and the 
parallel plate region can sustain at least one mode consistent with the polarization 
which we have employed. In order that a second mode not propagate in this parallel 
plate region, we must further assume that a <A. We also assume that there is a single 
reflected wave. Such a restriction puts further limitations on a/\ as well as on 8. 
These restrictions will appear when we have obtained the solution of the problem. 

3. Fourier transform solution of the integral equation. Before we turn to the 
Fourier transform solution of the integral equation (2.4) we shall first convert it into 
one of the Wiener-Hopf type. We note that the surface current density of the mth 
plate has the same magnitude as that of the zeroth plate provided we measure the 
distance along the mth plate from its edge. Hence, the surface current density on the 
mth plate differs from that of the zeroth plate only by a phase factor. This phase fac- 
tor arises because the amplitude of the incident wave differs from plate edge to plate 
edge by the factor 


eik(b cos §6+a sin 9). 
Thus 
I n(2 sitie mb) = To(z)e**#™¢> cos 6+a sin 6). 


where J(z) is the surface current density on the zeroth plate. Equation (2.4) may 


then be rewritten as 





: 28 Ks —- 
0 = dine(z, ma) + = pk f To(z)e‘*™H [Ry/(z — 2’ — mb)? + (n — m)*a?|dz’, (3.1) 
m=— 20 0 


where p=d cos 6 +a sin 0. If we replace z by 2 +b, Eq. (3.1) will read 


0 = etk ((2etnb) cos 6+na sin 6] 





1 - si ° " \ ae ee , 
+ r >» f To(z')et*omH [ka/} (m — m)b+(z— 2’) }? + (n — m)?a?|dz',z > 0. 
mM=—20 0 
Finally, when we divide the last equation by e***" and put m—n=gq, we get 


ib ve 7 . 
QO = eikzcos 6 4 : yi f To(s')e**°tH ) [ky/q?a? + (qb +2 — 2’)? |dz’ (3.2) 
q=—@ 0 





and this equation is of the Wiener-Hopf type. 
In order to put this equation into a form which amenable to solution by Fourier 


transform methods, we extend it for negative z to be 





¢:(z) = mM f To(2')e**#*H | ka/q?a? + (qb +2—2')?|dz’, 2 <0, (3.3) 
) 0 


i 
tion 


where ¢;(z) is an unknown function which is, save for a phase factor, the tangential 
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component of the scattered electric field at x=ma. In view of the periodic nature of 
the structure, the dépendence of the integral equation on m is not explicit. We may 
now replace Eggs. (3.2) and (3.3) by the equation 





o1(z) = do(z) + = 2 f To(2')e**#*H) [ka/q?a? + (qb + 2 — 2’)? ]dz’, (3.4) 


q=—2 


where now 
o,(z) =0 for z>0, 
Io(z) =0 for 2z< 0, 


0 for z <0, 
we = cos@ for z> 0. 

For analytical convenience, it is now assumed that k has a small positive imaginary 
part. This is tantamount to assuming that the medium is slightly absorbing. 

Before we can apply the Fourier transform in the complex plane to the solution 
of Eq. (3.4) it is necessary to study the growth order of the functions ¢,(z), Io(z) 
and ¢o(z). It is clear from a direct study of the integral Eqs. (3.2) and (3.3) that these 
functions are integrable for all finite z. The half planes of regularity of the Fourier 
transforms of $o(z), @:1(z) and Jo(z) are, of course, determined from their growth orders 
at infinity and we now proceed to determine these orders. Since we know ¢o(z) ex- 
plicitly, it is clear that its Fourier transform is 


ni 1 
—iwz z dz’ = 
J, — i[w — k cos 6] 


and is regular in a lower half of the w plane defined by the inequality {mw < Ym(k cos 6). 
Save for a translation on the z variable and a phase factor which is independent of z, 
Io(z) is, in certain units, the surface current density on any metallic plate. For z suffi- 
ciently large and positive, I9(z) is asymptotic to the surface current density in any of 
the parallel plate regions, that is, it is asymptotic to e**. Since Jo(z) is integrable at the 





origin, the Fourier transform of Jo(z), that is 


o 
f Io(2’)e~*”*'dz’, 
0 


is regular in some half plane defined by 


RekImk 
Qmw < Ym(K) ~ — | 7 
K 





> Imk, 


since Re(k)/| «| >1. 
We now investigate the asymptotic form of ¢:(z) for z large and negative. Before 
doing this, however, it is convenient to give another representation of the kernel of the 


integral equation (3.4). The kernel 





‘a y ea) | ky/q?a? + (qb + 2)?] 


q=-—-@ 
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has the Fourier integral representation 


© pikeatilalaY k*—w*—iwgd 

es eivz ) - dw, (3.5) 
4a Cc q=—oo a/ k? os w? 

where C is a contour which lies in the strip of regularity of the sum in (3.5). It is 
closed in the upper or lower half planes by a large semi-circle which passes between the 
poles of this sum depending upon whether z>0 or z<0. The strip of regularity is, 
of course, determined by the region in which the infinite series in (3.5) converges. A 
direct study of this series will reveal that the ordinates of convergence are given by 
the inequality, $mk cos (2a—86) <Ymw < Ymk cos 6. This now clarifies the reason why 
we imposed a small but positive imaginary part on k. Had we not done this, the series 
would only converge on the real axis of the w plane and as we shall see in the actual 
solution of the Wiener-Hopf equation, this situation would have presented us with 
some analytical difficulties. 

We may now write the sum in the integral (3.5) in closed form as 


i 








i f e sin an/k? — w? dw 
4nd co /k* — w[cos a\/k? — w? — cos (kp — wh) | 


For z<0, we close the path C in the lower half of the w plane. The poles in the lower 
half plane are w=k cos (2a—6) and two infinite sequences of poles both of which have 
negative imaginary parts. We shall have more to say about this double set of poles 
presently. Suffice it to be noted at this point, that the kernel has a second representa- 
tion which for <0 may now be written as 


eikz cos (2a—6) 





: + terms which attenuate exponentially for z large and negative. 
2ak sin (a — 6) 


It is clear then, that for z large and negative, $;(z) is asymptotic to 


f eik(2—-2") cos (2a—é@) 
Io(2’)dz’, 
0 2aksin (a — @) 





and thus, the Fourier transform of ¢;(z), i.e., 


0 
f e~ #9, (z)dz, 


—o 


is regular in the upper half of the w plane §mw> Smk cos (2a—6). 

The Fourier transforms involved in this problem then have a common strip of 
regularity, ¥m((& cos (@—2a)) <¥mw < Sm(k cos 6) and it is thus permissible to apply 
the Fourier transform to the integral equation (3.4) within this strip. 

Let ®,(w) be the Fourier transform of ¢;(z) and J(w) the Fourier transform of 
' Io(z). The Fourier transform of the integral equation (3.4) is then 


1 J(w) sin a/k? — w? 
®,(w) = + = —— —- =: (3.6) 
i(w—kcos0) 2/k? — w* [cos ar/k? — w? — cos (kp — wb) ] 





The Wiener-Hopf theory now tells us that we can split this transform equation into 
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two parts. One part will be regular in an upper half plane, ¥mw> mk cos (@ — 2a), the 
other in a lower half plane §mw<Qmk cos 6 and both of these half planes have a 
common region of regularity. It is well to note here that we use the term regularity 
in a slightly extended sense. We imply by regularity that the function has neither 
zeros, branch points nor poles in the region of regularity. That is, the function as 
well as its reciprocal is “regular” in the conventional sense of the term. Suppose we 
assume that we can write 


K_(w) sin a\/k? — w? 
K,(w) i o/k? — w? [cos ar\/k? — w? — cos (kp — wb) ] 





where K_(w) is regular in the proper lower half plane and K,(w) is regular in the 
proper upper half plane and that there is a common strip of regularity for both K_(w) 
and K,(w). Then 





J ” Ki(w) | J(w) K_(w) 
teed i(w — k cos 8) 2 te 


The left side of Eq. (3.7) is regular in an upper half plane while the second term on 
the right side is regular in a lower half plane. The term 


K,(w) 


i(w — k cos 6) 





is only regular in the strip of regularity. This function may be decomposed into two 
functions in such a manner that one function is regular in the appropriate upper and 
the other in the appropriate lower half plane, since 


K(w) 7 K,(w) — K.(k cos 6) K(k cos @) 


i(w — k cos 8) 5 i(w — k cos @) i(w — k cos @) , 





The first term on the right no longer has a singularity at w=k cos 0, but is regular 
in the upper half plane and the second term is regular in the lower half plane. Thus 
Eq. (3.7) can be rewritten in the form 


K.(w) — K,(R cos @) _ J(w) K_(w) K(k cos @) 


» (3.8) 
i(w — k cos 8) 2 i(w — k cos @) 








#;(w)K,(w) — 


The right side of the equation is regular in the lower half plane ¥mw< mk cos 6 
while the left side is regular in the upper half plane §mw> mk cos (@—2a). Both 
sides have a common strip of regularity and hence the left side of (3.8) is the analyti- 
cal continuation of the right side. Such an equality can only hold if both sides of 
Eq. (3.8) are equal to an integral function, that is, a function regular everywhere in 
the complex w plane. We have then 


J(w)K_(w) K(k cos 6) 


. = integral function (3.9) 
2 i(w — k cos @) 





and also 
K,(w) — Ki(kcos@) ; 
©,(w)K(w) — = integral] function. (3.10) 
i(w — k cos 6) 
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We shall now show that it is possible to decompose the function 


sin ax/k? — w 


/ ee — w v2 [cos ax/k? — w? — cos (kp - _ wb) | 
into two functions, one of which is regular in the lower half plane §mw<Qmk cos 6, 
while the other is regular in the upper half plane {mw>Qmk cos (@—2a). The de- 
nominator of the fraction may be written as 


cos a\/k? — w? — cos (kp — wh) 
ee wth-wl . &—-o~ws/Fo 


sin - sin 
. 2 


= [a/R — w+ kp — wh] [kp — wh — ar/P? — w*] 





ll 
bt 











~ (ar/k? — w? + kp — wh kp — wh — a {ali i 
x TI{1-** Ba “) [1 - (kp an ) 
anare 4n?r? ae 4n?r? 
ay/k? — w? + kp — wh — 
= 3[(kp — wb)? — ab? eT anes on eto teeta 
n=1 -— 





x TI 


n=—oo 


x Tl 
x I 


n=—eo 


Sa at hy —_ 
}1- “. ih om + p =e (a™ k?—w*+ kp—wb) /2nx 


one 





kp — wh — ar/k? — w? —_— 
E — . V Jew ~wb—a k?—w?) [2ne 
2nr 
kp — wh —ar/k?— w —, 
}1- = alte =) e(ke—wb—aY k?—w?) /2ne 
2ne 
The exponential factors in each of these products has been inserted to render the 


products absolutely convergent. The above expression may now be rewritten to read 


. mon’ kp — wh)? a%(k*? — w’) 
$(a? + b*)(w — o1)(w — o2) TI it = +=) “ 22 ew eater; 4.80) 
4n*x? 


eas 2nr 





where the prime on the products denotes the absence of the term m =0 in the product. 
The infinite product in the last expression may now be expressed in a manner such 
that it puts into evidence the portion which is regular in the correct upper half and 
lower half planes. Indeed we may express (3.11) as 

i iV, Jel (tow b+ wai) [240] +i(4/2-a) 


1(a? + b*)(w — o1)(w — o2) []’ [A, 


n=—o 


ce) 
x Il’ [A, + £¥,, el (te b—-wast) /209}—i(2/3-2), 


n=—0oo 


where now 
o1 = kcos8, oo = k cos (2a — 8), 


and 


kp \? ak \? wa CSC @ 
A, = 4/ sin a (1 ba ) _ ( ) ’ Y,, = cosa 1 a ) + Se Te 
2nn 2rn 2an 2rn 
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where again, the exponential factors following the infinite products have been chosen 
to insure the absolute convergence of the product. One should note at this point that 
the choice of these exponential factors is not unique and indeed need only be asymp- 
totic to the factors which we have chosen. However, we shall see that a second in- 
tegral function x(w), introduced into the decomposition of K(w), is determined in 
terms of the factors which we have chosen. We have finally that the factor 


(w seit 71) TI [A,, all iW, Jette obvi) 2enticelta) T] [A, a iV, Jel te—wb—wai) /24n]— i(r/2—a) 


—o n=1 


has no zeros in the lower half plane ¥mw< mk cos 6, while the factor 
p 
oo —1 
(w — ge) II [A,— iW, Jel -w b—wai) /24rn]+i(4/2—a) Il [A, te iV, Jel te—wb-wat) /2en]—i(x/2-a) 


n—l n=—o 
has no zeros in the upper half plane §mw>9m[k cos (@—2a) |. The factorization of 
sin a/k? — w? 
/ k? — w 


is more direct, for 


sin a\/k? — w? a*(k? — w? 
het oT[1~- ait | 
a/ k? — w? ani n*x? 
a’ ; . iaw] 
ett (w “a jew tev! (wy oo x)eiaw!* ] 1— (= 1-(£) oe oa e7 iawirn 
T n=2 rn wn 
=. ~ fak\? iaw 
* | 4/1 - (< -) au = etal rn. 
n=2 wn 


a = , _ (2k? 3 Pa 
eens (w che k)e -iew/* T] e7iaw/rn 
us n=2 lL awn 


has no zeros in the lower half plane Jmw<Jmx, sada ha factor 


a aie ak\? iaw] | 
race (w f xjeiaw/* T] 1 a Ree ees egiaw/an 
v n=2 7 Nn 


has no zeros in the upper half plane Jmw>Qm(—x«). We thus find that 


1aw 
[1 —_ 1-(4)4 +] remem = (w om x)e~taw!™ eX (w) 
en Nod 
K_(w) = — - 


(w — o1) Il [An se iVp)el(to-wbtwai)/2en}ticr/2-a) T] [An + EW,, Jel? -wb—wai) /2xn}—i(*/2-2) 


N=—0o nal 


is free of zeros and poles in the lower half plane $mw< {mk cos 6. The factor ex™ 
will be determined so as to make K_(w) have algebraic growth as |w|—>© for 
¥mw <0. With x(w) so chosen, the integral function sought can only be of algebraic 
growth for |w|—«. K_(w) is regular in the lower half plane ¥mw< mk cos 0. 
Finally, 

(a? + b*)(w — ex] ] [An — EV, Jel eo -whtwas)/2n}+i(x/2-a) Il [An + iW, Jel 4? -wb+ wai) /2en}—i(*/2-) 





II 


The factor 
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K,(w) = - anne 
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2a Il [/1-(2)- — — — | giaw/tn ll + x)eiaul™ 
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has no zeros or poles in the upper half plane ¥mw> mk cos (@—2a). 

We shall now discuss the asymptotic form of K_(w) as | w| —o, ¥mw<0. This 
procedure will enable us to determine the unknown integral function x(w). It has 
been shown by Schwinger’ that functions of the form of K_(w) are independent of ka 








for |w|—-o, ¥mw<0, and r<ak<2zx. Thus 
eX wv) II [1 ao aah e7iaw/an bed e tau/* 
na? mn T 
7 a 





2an 


= ; wa CSC a 

Il sina — ¢[ cos @ + ——— } | el(wot-we)/8en]+i(4/2-2) 
or saemeeaie (3.12) 
)] e7 [(waitwb [/24n]—i(*/2-a) 








"| ; 4 wa CSC a 
II sin a + i{ cos a + 


n=l 2an 


The products in (3.12) are now in the form of gamma functions and 


wa csc a\? — wa csc ae~** wa csc ae‘ 
ex (w) eivayl« Tr r panacea 
2r 2 2a 
iaw iaw\ - iaw 
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T T Tv 


where ¥ is the Euler-Mascheroni constant. Using the Stirling expansion theorem 
—o, ¥mw<0 we get 








K_(w) ~ 


P 








for 











WwW 
aw CSC a —[(wa exe @)/2™Je—**—1/2 7 gw csc ‘ 
eX wg csc? a (- a “1, ci) -— ) el(wa csc a) /2}e**—1/2 
T T 
K_(w)~— —_— ila 
: iaw\ “aw!*)-1/2 
4% (=) 
TT 
Cex wv) tiaw/* [(a-*/2) cot @+In (csc &)/2) 
~~ ’ 
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where C is a constant. Thus if we choose 


— taw 7 
x(w) = | (« - =) cot a — In2sin «|, 
T 2 


K_(w) will have algebraic growth for | w| large, ¥mw <0. 

Now J(w), which is proportional to the Fourier transform of the surface current 
density on the various plates, approaches zero for | w| large, §mw<0. This assumes, 
of course, that J(z) can at most be of exponential growth for 2 large and positive and 
is integrable for z finite. Thus K_(w) J(w) approaches zero for | zw | large and ¥mw <0. 
If we now return to Eq. (3.12) we see that as | w| becomes large, ¥mw <0, the integral 
function in (3.9) is asymptotic to zero. We may now apply the same argument to Eq. 
(3.10) and find that the integral function is again asymptotic to zero. But by a theo- 
rem of Liouville, and analytic function which is bounded in the entire complex plane 
is constant and in this case the constant must be zero. We thus have 


21K(k cos 6) 
 K_(w)(w — k cos 8) 
If we were interested in the explicit form of the surface current density, we could 
obtain it from J(w) by evaluating the Fourier inversion integral 








J(w) 
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~ f K(k cos 0)e™*dw 
aJc K_(w)(w — k cos @) ' 


where C is a contour which may be taken as a straight line within the strip of regu- 
larity of the Fourier transforms of I(z), ¢:(z), @o() and K(z). The contour is closed 
above by a semi-circle, which by familiar arguments in contour integration may be 
shown to make no contribution to the value of the integral. In the next section we 
shall show that it is possible to find the reflection and transmission coefficients with- 
out evaluating this integral in detail. 

4. Investigation of the far fields. In order to find the reflection and transmission 
coefficients, we now investigate the asymptotic form of ¢(x, z) for |z| large. To this 
end we note that Eq. (2.3) can be written in Fourier integral representation as 





i bed eieztik mp—iw mb+|2— mal k*—w* J(w)dw 
:, 2) = Pine(x, 2) + — 7 
(3 ) e ( ) t 4nJec 27 Vk? —— 


where C is the contour which we described at the end of Section 3. This in turn, may 
simplified to 
o(x, 2) = dine(X, 5) 

[sin Vk? — w? (x — an — a) + ef (4?) sin \/k? — w? (an — x|dw 


i 
a ag) \piwati(kP—wb) - a 4.1 
re J (whe Vk? — w* [cos a\/k* — w? — cos (kp — wh) ] Bay 





where m is the greatest integer contained in x/a. From (4.1) one can get the asymp- 
totic form of }(x, 2) as z becomes large and positive. Since J(w) is regular in the lower 
half of the w plane {mw< mk cos 6, we can close the contour C by a large semi- 
circle which passes between the poles in the upper half plane. For na<x<(n+1)a it 
can be seen that due to the form of the integrand, there is no contribution from this cir- 
cular arc as its radius becomes infinite. In the upper half plane ¥mw > ¥mk cos (2a—8), 
there are two poles which correspond to propagating modes, namely w=k cos @ and 
w =x. All other modes are.attenuated modes in the sense that they have large positive 
imaginary parts compared to the imaginary parts of k cos 6 and x. If we now express 
J(w) as a function of w and use the above described contour in the evaluation of the 
integral in (4.1) we have then to consider the asymptotic form of 

s f cilke—wo]ntiws [sin (x — an — a) Vk? — wt + eft? sin Vk? — wt (an — 2) ]Ky(k cos ew 
2rlc (w — k cos 0=)K,(w) sin aV/k? — w? 

This in turn is equal to [e**( sin ¢+2 cos &) _ Tei sin xx/a +terms which approach 
zero for 2>>0]. For z large and positive, this is asymptotic to 





Wx 
Pine(x, 2) — Te*** sin — - 
a 


Hence, save for a numerical factor, the functional form of ¢(x, 2) as z becomes infinite 
is e*** sin rx/a, that is, it represents a travelling wave in the parallel plate region with 
propagation constant x, as it should. The amplitude of this wave is 


mein(ko—«)(—) [1 +4 ei(te—«b)] K(k cos 8) 
(x — k cos 0)a%K4+(k) 





T=|T|e® = 
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and depends of course on the particular parallel plate region for which it has been 
computed. Since T is the amplitude of the wave transmitted in the parallel plate re- 
gion it is the transmission coefficient because the amplitude of the incident wave has 
been taken to be unity. If we now assume that & is real, the magnitude of T is 


23/2 sin (a — @) 


= oanseeaapeipan =f 


~ -4/(k cos 6 + «)(« — k cos (2a — 8) 


a 








a quantity independent of the particular parallel plate region considered. Its phase 
angle depends, of course, on the particular parallel plate region. We shall not give 
the phase angle explicitly since we shall not use it in our later discussions. 

For z large and negative we close the contour in the lower half of the w plane. 
There is again no contribution from the circular arc which is drawn between the poles 
in the lower half plane and so we need only evaluate the residues from the poles in the 
lower half plane. The dominant contribution now arises from the pole w =k cos (@ — 2a) 
and in this case the dominant term is 


K(k cos 6) et* lz sin (2a@—@)+z cos (2a—@)] 


k[cos (2a — 6) — cos a] K, [k cos (2a — @)] 





all other terms in the integrand approaching zero for z large and negative. Here 
K,! [k cos (2a—6) | means, as usual, the derivative of K,(w) with respect to w evalu- 
ated at w=k cos (2a—6). The amplitude of the reflected plane wave is the reflection 
coefficient R if the amplitude of the incident wave is taken as unity, so that we now 
have 


K(k cos 8) 
k[cos (2a — 0) — cos 6|K,’ [Rk cos (2a — 6) | 





R= 


Assuming, once again that & is real, the reflection coefficient may then be rewritten 
in complex polar form as follows: 





R too 4 / (2a — 0) + x) 
= e* 1 2 ] 
(k cos 6 + x)(k cos (2a — 6) — x) 





where now 








( ke ] 
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, ka cos 6 ka 
arc sin . —_——________— — — cos 9 


ka\? | ™ 
Nn 1—(|{—)} sin*é 
t 7m 
ak ri 
— — cos oft + (« — =) cos a — In 2 sin af : 
T 2 


Ms 


+ 


[<) 


n= 





and 
p | cosa + sin (0 — a) | 
; 2an ka T 
@.= — >> jarcsin — Serene — — cos (@ — 2a) —{— — «) > 
oat ka . 2nr 2 
1 — — sin (2a — 6) 
™m J 
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} cos a@ a — sin (6 _ a) 
re > ; 2an ka ( 2a) (< ’ 
‘arc sin - — — —— cos (0 — -({—- 
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| i 1- — sin (2a — 8) | 
7 
( ka 
—cos (2a — 6) 
ae A , 7m ‘ ak 
+> {arc sin —_____________. — — (cos (2a — 6) 


mee | / ka\? . nT 
! — ~) sin? (2a — 6) 
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ak { T . 1 
— cos (2a ~ 6) 41+ (a-) cosa ~ In 2sina ; 
© 2 f 

It is evident that there will be restrictions on a/d and 6 if we are to have a single 
reflected plane wave. These restrictions become evident when we study the arc sin 
sums and observe that conceivably the first term in the sums beginning with index 
unity can exceed unity. We have tacitly assumed that they do not, for otherwise 
they would appear in the amplitude factor as real terms. Thus we must see what is 
implied by the condition that all factors in the infinite products be complex, or equiva- 
lently A2>0. If we demand that Aj>0 it is clear that all other A,s, n=1, 2, - -- will 
also be >0. The condition A;>0 is equivalent to 








ak 2a sin @ 
(i) —=—< 
7 r cos? $(@ — a) 
and 
. ak sin a 
(11) — _ 
7 sn? $(@ — a) 


Condition (ii) is always satisfied since a/ is always positive. Condition (i) can be 
more restrictive than the condition 1/2 <a/A <1. For example, if @=2/12, a=52/12, 


then condition (i) implies 
a/\ < .65. 
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5. No propagation in the parallel plate regions. In the Fourier transform solution 
of the integral equation (3.4) we have assumed that there was only one propagating 


mode in the parallel plate region, i.e., 
_ , w2 
I n(z) ~ e*** sin —» k>0O 
a 


for z large and positive and z in the parallel plate region. Suppose now we dispense 
with this assumption and ask what form the reflection coefficient takes if we now as- 
sume that x?<0, i.e., 0<a/A<1/2. In this case x is, of course, imaginary and 

wx 


In (2) ~ EY 1a) —# ¢ sin — 
a 


for z large and positive and z in the parallel plate region. The result we desire can be 
obtained most easily by studying the result which we have obtained in Section 3. 


We note that if « is purely imaginary and k is real, the amplitude of the reflection 
coefficient becomes complex of magnitude unity. Indeed for x?<0 





(k cos @ — x)(k cos (2a — 0) + x) 
(k cos 6 + x«)(k cos (2a — 6) — x) 





f/f 


ka cos 6 4 ka cos (2a ~~ 9) 
— arc sin : 


ef “ (2) sito sin? 0 a oe 1— (*)'sint sin® (2a — 6) ( 


Thus for this situation, the amplitude of the reflection coefficient is —1. The phase 


angle @j is given by 











sad om} arc sin — 








ka cos 8 
e/ = 0, + arc sin eae 
ka\? . 
“4/1 — (~) sin? 6 
us 
while @/ is now given by 
_ ka cos (2a — @) 
@/ = 0. + arc sin — —,, — 


T ak 
A — (=) sin? (2a - - 6) 


Hence the reflection coefficient for 0<a/A<1/2 is now —e*@1’-®’), For a single re- 
flected wave, the inequality (i) in Section 4 must still be satisfied, although now it is 
not as severe. 

6. Discussion of results. It should be pointed out that some of the results ob- 
ained from our calculations can be interpreted in a simple physical manner. For con- 
venience, in this discussion, instead of the angle 0 we use the angle 7, which the inci- 
dent wave makes with the normal to the trace of the edges of the plates. It is readily 


verified that 
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and also that the angle r which the reflected wave makes with the normal is also 
equal to z. The condition that there be only one reflected wave 


2a sin a 


d cos? 4(@ — a) 
is seen to be a result of simple grating theory. If the waves scattered by a uniform 
grating are not to interfere constructively in the region from which the waves are 
incident (except for the specular case r=i) the condition a’(1 +sin 7)<A must be 
satisfied where a’ is the distance between neighboring scatterers. In our case 
a’ =a csc a. Expressing @ in terms of 1, the relation 





2a sin a 


r cos? $(@ — a) 





is seen to be equivalent to a csc a(1 +sin ¢)<X. If this condition is satisfied and the 
condition for no propagation, \> 2a, is also satisfied, the plates act as a perfect plane 
mirror. However, while the magnitude of the reflected wave is unity, its phase is not 
mt but O{ —@/. It is easily shown that it will be 7 on any plane parallel to that of the 
trace at a distance d given by 

(4rd/X) cos i + 2mr = OF — Of m=0,+1,+2,---. 


Therefore, as far as all distant fields are concerned, the plates behave in this case like 
a perfect plane mirror whose surface coincides with any of the planes given by the 


above equation. 
When transmission is possible in the parallel plate region the wavelength in this 


region differs from that in free space. One would, therefore, expect to find some anal- 
ogy with the phenomena associated with a plane interface between two dielectric 
media. This can be shown for the case a=7/2. In this case the magnitude of the re- 
flection coefficient is 

kcosi—k 

| Be a cee 

kcoosi+x 
This expression is identical with that obtained for the reflection at a dielectric inter- 
face of a wave with the electric vector parallel to the interface. The phases are differ- 
ent in the two cases and one can again find a set of planes at a distance d from the 
trace given by 

(4rd/d) cos i + 2mm = O2 — Oy, m=0, + 1,--- 


such that the distant fields are identical in the two cases if we regard any one of the 
planes as the interface. 

The expression for the magnitude of the reflection coefficient should be of use in 
estimating the reflection of waves incident on a metal lens provided that the radius 
of curvature of the lens (i.e., the angle a) does not vary too rapidly. 
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A PROBLEM IN THE PROPAGATION OF SHOCK* 


BY 


MONROE H. MARTIN 
University of Maryland 


Introduction. This paper deals with a single problem in the rectilinear motion of 
a gas, namely, what is the subsequent behavior of a gas initially at rest if its initial density 
is a constant po in the region |x| <1 and a constant p2<po in the region |x| >1? 

The behavior of the gas is an idealization of the behavior of the atmosphere in 
an infinitely long right circular cylinder after an explosion within the cylinder. 

It is assumed that the pressure p and density p of the gas are related by the isen- 
tropic law p~=k?p7 where k? is constant for all x and all ¢. Under the law of conserva- 
tion of energy! (Rankine-Hugoniot equation) there is a change? in entropy across a 
shock and the results in the paper may be regarded as an approximation to the actual 
state of affairs only in the case where the change in density across the shock is very 
small with a correspondingly small change in entropy. 

At times the author has not hesitated to restrict attention to a monatomic gas 
(y =5/3) in order to avoid formal mathematical difficulties.* The behavior of the gas 
undergoes marked changes as the difference po — p2 between the initial densities is per- 
mitted to vary.‘ 

If po— 2 is sufficiently small the two initial shocks give rise to shocks traveling in 
opposite directions towards infinity as ¢ increases indefinitely. Up to a certain instant 
the shocks travel with constant velocity greater than the velocity of sound in the un- 
disturbed gas. After this instant their velocity of propagation decreases monotonically 
with time, to approach the velocity of sound in the undisturbed gas as the shocks 
recede to infinity.’ The behavior of the gas between the two shocks is followed up to a 
stage when the mapping® of Riemann’s (r, s)-plane upon the (x, ¢)-plane loses its one- 
to-one character. The further behavior of the gas still awaits determination. 

Plates 1 and 2 at the end of the paper present qualitatively the variation of density 
(or pressure), over the gas for po— pe sufficiently small. 

1. Fundamental principles. Assuming that the pressure p is a monotonic increas- 
ing function of the density p and denoting the velocity by u, the partial differential 


equations 


* Received March 8, 1946. 

1See Riemann-Weber, Die partiellen Differentialgleichungen der Mathematischen Physik, 6th ed., 
Friedr. Vieweg & Sohn Brawnschweig, 1919, vol. 2, pp. 549-550. 

2 Indeed under the Rankine-Hugoniot hypothesis it follows from formula (10) on p. 513 of Riemann- 
Weber, op. cit. that the entropy of the gas in back of the shock depends upon the ratio of the densities of 
the gas on the two sides of the shock. This ratio changes as the shock propagates and consequently the 
the entropy of the gas in back of the shock is not constant. 

8 The examination of other values of 7 has been begun by R. C. Rand in his doctorate thesis entitled 
The rectilinear motion of a gas subsequent to an internal explosion. A copy of this thesis is on file in the li- 
brary of the University of Maryland. 

*R. C. Rand, loc. cit. 

5 See, however, the last sentence in §6 of the present paper. 

6 For a discussion of this mapping see Riemann-Weber, loc. cit., pp. 533-536 or The rectilinear motion 
of a gas, Amer. J. Math. 65, 391-401 (1943). This paper will be cited as I. 
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p(u,+ ust) +G%p.=C, prt (pu), = 0, G? =G%) = 7’, 
for u, p become 
r,t ar, = 0, si + Bs, = 0, a=u+G, B=u-G, (1) 


if we set 
r) 
u=r+s, f (G/p)dp =r—s>O0O. (2) 
0 


Clearly u, p are monotonic increasing functions of r+s, r—s respectively and 
a=al(r,s)=r+s4+G(p(r—s)), B= Br, 5s) =r +s —G(p(r — 5)) (3) 
satisfy 
a(— s, — r) = — A(r, 5), B(— s, — r) = — a(r, s). (4) 

A point of the (u, p)-plane, or its correspondent by (2) in the (r, s)-plane, is said 
to represent or be a state of the gas. The points of the (r, s)-plane representing states 
of the gas comprise a half-plane r2=s termed the state plane. Points representing states 
having the same velocity (density) lie on the lines r+s=const. (r—s=const.) and 
the velocity (density) of a state increases with the distance of the point (r, s) from 
the line of zero velocity r= —s (the line of zero density r=s). The velocity is positive 
or negative according as (r, s) lies above or below the line r= —s. 

In general a solution r=r(x, t), s=s(x, t) of (1) transforms a region of the (x, ¢)- 
plane into a region in the state plane and is single valued. The inverse transformation 
T: x=x(r, s), t=t(r, s) is not necessarily single-valued and is regarded as assigning 
the state (r, s) to its transform (x, t). Corresponding to (1) there is the system 

x, — Bt, = 0, x, — al, = 0, (S) 
of partial differential equations for x(r, s), t(r, s) in T. The Jacobian J of T is 
J = — (a — B)td, = — 2Git,A,. (6) 
If x, ¢ are solutions of (5), the system of Pfaff 
x—at—v x— Bpit-—o 
dw = (x — at)dr + (x — Bé)ds, dv = 2 ——_————_ dr — 2 ds, 
a-B a-—-8B 
is completely integrable, and conversely. When we write 
x — at = w,, x — Bi = w,, (7) 
the integrability condition for the second equation becomes’ 
(a — B) Wr. — B,(w, — w,.) = 0. (8) 
Taking w a solution of (8) it follows from (7) that a transformation T is 
Bw, — aw, Wr — We 
Ym s= — ————_) t= — ———__-- 
a-—f8 a-B 
The following theorem is a direct consequence of (4). 


7 Cf. Riemann-Weber, loc. cit., pp. 536-538 or pp. 393-394 of I. 
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THEOREM 1. Given w=w(r, s) a solution of (8), another solution is ®=w(—s, —r) 
and T., Ts map points which are reflections of each other in the line of zero velocity 
r=—s, into points which are reflections of each other in the line x =0. 


As a corollary, we see that if w(r, s)=w(—s, —r) points which are reflections of 
each other in the line r = —s, are carried by 7, into points which are reflections of each 
other in x =0. 

The theorem is obvious a priori on physical grounds. Given any motion of the 
gas, its particles may be reflected in the plane x =0 to gain another motion. 

Taking r=ro=const., the second equation in (1) upon multiplication by d8/ds 
becomes 

a(x — Bt, B) _ 


B: + BB. = 0, 
0(x, t) 








and therefore a solution of (1) is given implicitly by® 
rm te, a Bt = v(8), (9) 


W(8) denoting an arbitrary function of 8. Corresponding to s =s9=const., a solution 
of (1) is obtained from 
x—at= (a), S$ = So. (10) 


For a fixed s in (9) the state (ro, s) is assigned® to all points of the straight line 
x —Bt=W(B8). This line is termed a propagation line and the state (ro, s) is said to be 
propagated along it. Physically the state (ro, s) is propagated through the gas with a 
velocity 8 with respect to a fixed plane. 

Let us assume that 7, puts the states of a region R of the state plane in (one-to- 
one correspondence with the points of a region X of the (x, ¢)-plane. The transform 
by 7. of a segment of r=const. (s=const.) in R is a curve in X termed an r-curve 
(s-curve). The r and s-curves provide a curvilinear coordinate system on X from which 
the state of the gas may be read off at any point of X. 

From (5) the slope of an r(s)-curve’® is 1/a (1/8); from (9), (10) the propagation 
lines drawn from the points of an r(s)-curve have slope 1/8 (1/a). Therefore the 
tangents drawn to s(r)-curves at the points of an r(s)-curve are propagation lines and, in 
so far as they do not intersect, may be used to assigned the states on the r(s)-curve to the 
points of the region covered by them. 

Two r(s)-curves C, C transforms of r=ro (s=5So) under Ty, Ts respectively are 
said to be propagated from each other if the propagation lines drawn from points of 
C, C which are transforms of the same state are identical. 


Lemma 1. Two r[s]-curves C, C transforms by Tw, Tz of r=ro [s=so| are propagated 
from each other if, and only if w,(ro, s) =W.(To, 5) [w,(r, So) =@,(r, So) |. 
From (7) parametric equations of C, C are 
C: x — a(ro, s)t = w,(ro, 5), x — B(ro, s)t = w,(ro, 5), 
C: x — a(ro, s)t = W,(Tro, 5), x — B(ro, s)t = (To, S). 


8 Cf. Riemann-Weber, loc. cit., p. 518. 
® Cf. Riemann-Weber, loc. cit., pp. 516-520. 
10 First noted by R. C. Rand. 














1947] A PROBLEM IN THE PROPAGATION OF SHOCK 333 


Along a propagation line propagating the state (ro, s) we have x—{(ro, s)t =const. 
Hence propagation lines drawn from a point on C and a point on C, both transforms 
of the same state (ro, s) will be identical if, and only if, w.(ro, s) =W,(ro, $). 

It is interesting to note that tangents drawn to C, C at points which are trans- 
forms of the same state are parallel. 


LemMA 2. Given two r[s]-curves C, C which are propagated from each other, curve © 
will pass through a point (&, 2) on the propagation line propagating the state (ro, 5), 
[(¥, so) ] from C if and only if @,(ro, 3) =£—a(ro, 5)i[w.(#, So) =¥—B(F, So)é]. This con- 
dition determines w,(ro, s) [w.(r, so) | uniquely. 

The first part of the lemma follows from the parametric equations of C. To prove 
the second part we set r=7o, w=@ in (8) to obtain an ordinary, linear differential 
equation for #, since #,=w, is a known function of s. This determines #, uniquely, 
for w, is known when s =5. 

2. Shocks and buffer waves. Under the assumption that G increases with p it 
follows that G=G(p(r—s)) is an increasing [decreasing] function of r[s] for fixed 
s|r]; from (3), one concludes that a[] is an increasing function of r[s] for fixed s[r]. 


Lemna 3. If initially r=ro for —~ <x<+ and s=s, or Ss. as x<0 or x>0 with 
5; <So, subsequently the state of the gas is unchanged exterior to the “buffer region” between 
the lines x =B(ro, s1)t, x =B(ro, S2)t. Within this region the state (ro, s) with s,<s <se ts 
propagated" along the propagation line x =B(ro, s)t. 


Initial states are propagated along the propagation lines 
x — B(ro, :)t = ky < 0, x + B(ro, S2)t = ke > 0, B(r0, $1) < B(ro, 52), 


which diverge as shown in Figure 1 to assign the state (ro, s;) to the region on the left 
of OA, and the state (ro, s2) to the region on the right of OA». To obtain the states in 


Az 


\\eb LL, 


ANN \ vi 


@) 
Fic. 1. 











the buffer region A,OA2 one sets ¥(8)=0 in (9) and draws the propagation line 
x =B(ro, s)t from 0. Along this propagation line the state is (79, s) and as s ranges from 
S; to Sp the propagation line turns from OA, to OA; to assign states to all points of the 


1 Cf, Riemann-Weber, loc. cit., pp. 520-521. 
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buffer region. It will be observed that the states vary continuously along a line 
t=to>0. 

In this solution of (1), the inverse transformation T is not single-valued the, seg- 
ment sis S52 of r=ro being carried by T into the half-plane t20. 

Physically the buffer region corresponds to a disturbance P,P» affecting two bodies 
of gas of different uniform states in contact with one another initially, the end points 
of the disturbance traveling with the local velocity of sound in the two bodies of gas. 
The passage of this disturbance through the gas is termed a buffer wave. 

A shock exists at x=& if p1%p2 and is propagated with a velocity’? 


lon = De lan Os — Oe 
jondy/S Slew gfe Sa ol (11) 
Pi Pi — pe Ps Pi" Pe 
where 
m= ulE—0), m=e(E—0), pr = p—E—O), 
uz = u(é + 0), p2 = p(§ + 0), pe = plé + 0). 


The curve x= &(t) in the (x, #)-plane is termed a shock curve. It will be sufficient 
for the purposes of this investigation to consider progressive condensation shocks 
arising when p:> 2 and the positive sign is taken in (11). For a shock of this type one 
has the condition 





My — U2 = V(pi — p2)(ps' — pi), (12) 
with 
— = (wip; — U2p2)(pi — po). (13) 


If (71, 5:1), (Y2, S2) denote the correspondents of (u1, p:), (ue, p2) by (2) and the state 
(re, Se) on the right of the shock is given, the state (71, s:) on the left of the shock is not 
uniquely determined but, by (12), may be any point of the curve. 





r+s=n+5.+V(p — poor! — op"), p= P(r—s), p=er—s)>p2 (14) 


in the state plane. This curve is termed the compatibility curve of the state (r2, 52) 
and its equation may be written in the parametric form 





r=Brntsstot V(p — pa)(ox' — pf, 


1f — o~) (14) 
s=3{r2+s2—v+V(p — po)(or? — p)}, = > 0, 





upon introducing the parameter v=r—s, where, of course, v2 =f2—Se. 


LemMMA 4. The compatibility curve of a state (r2, s2) rises with increasing r from the 
point (re, Se), at which it has a horizontal tangent. 


Both derivatives of 7, s with respect to v will be positive provided 


~ G%p) > 2 
p2 p — po 


= G*(p) where p2<5 <p, 





12 See, for example, Riemann-Weber, op. cit. p. 513. 
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which inequality holds for p>pe t 
under the assumption that G in- 
creases with p. 

We shall now consider the si- 
multaneous generation of a shock 
and buffer wave™ as pictured in B 
Figure 2 where B is a buffer region 
between two regions Ro, R; of uni- 
form state (7o, So), (f1, 51) respec- 
tively and OS is the shock line O 
separating R, from the region R» Fic. 2. 
of uniform state (re, Se). 








Lemna 5. A shock and buffer wave are generated simultaneously" at the contact of two 
bodies of gas of different uniform states (ro, So), (T2, S2) provided the point (ro, So) im the 
state plane lies directly underneath the compatibility curve of the state (re, s2). 


Choosing the state (r2, se) in Re arbitrarily, the state in R, must be represented by 
a point on the compatibility curve of the state (re, s2); and if this point lies directly 
above (ro, So) the existence of the buffer region B is assured by Lemma 3. 

3. The isentropic case. Here p=k’p? with k, y>1 constants and G increases with 
p so that the results of §2 remain in force. Moreover 





—] 1+ 3- 3- i+ 
ce ——(7 — $), (one aor Pee i (15) 
2 2 2 2 
and (8) becomes 
3-Y¥ 
(r — S)Wre — m(w, — w,) = O, m = ————— (16) 
2(y — 1) 


For monatomic gases y =5/3 one finds a=3(2r+s), B=3(r+2s). Also m=1 and 
(16) becomes 
(r — S)Wre — (wr — w,) = 0, (16’) 
the general solution of which is 
R-S 
w= , R= R(r), S = S(s), (17) 


y =F 





R(r), S(s) being arbitrary functions. The transformation T, is 








x= _ + 2s)w, — (27 + 5) Pa _3 wm as) 
es 2r-s 
or, from (17) 
2s) R’ 2 Ss! 
x= 3 ih (R o S) _ (r + s) + + 5) ; 
tien (r—s (18’) 


3(R — S) 3 R+S’ 


(r — s)8 2 (r — s)?” 


13 Cf, Riemann-Weber, loc. cit., pp. 527-529. 
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so that 


3 (r — s)?R” — Ur — s)(2R’ +S’) + 6(R — S) 








2 (r — s)4 
; (19) 
3 (r — s)2S” + 2(r — s)(R’ + 28’) — 6(R — S) 
os 2 ty a s)4 ] 


primes denoting differentiations of R, S with respect to their arguments. 


LemMA 6. For a monatomic gas the compatibility curve of a state (re, S2) is an arc of 
an algebraic curve of eighth degree ending at (re, s2), about which the points of the arc per- 
mit the expansions 


ee et ee ee a eee (ne twa es >. g = Ife. (20) 


At all other points r is a regular analytic function of s with a positive derivative and r, § 
are regular analytic functions (14') of the uniformizing parameter v, with respect to which 
they possess positive derivatives. 


For monatomic gases equations (14), (14’) become 





rT + 5S = 19 + Se + Vv Ps [(r on s)§ —_ (re “~ S2)*] [(re2 wae 52) * — (r — s)~4], (21) 
nt stot Ve = Der — oF), 


(ro + s2 — 0 + Vpg(v' — 0.) (or? — 0-4), 





f=: 


(21’) 





Nie tole 


5 = 


from which the statements in the lemma follow straight forwardly. 

We return to the general adiabatic case. A point in the state plane represents a 
state for which the velocity is subsonic, sonic or supersonic according as the point 
lies in, on the boundary of, or exterior to the region a>0, 8 <0 between the straight 
lines a=0, B=0. 


Lemma 7. The angle of inclination @ of the tangent at a point of an r{s|-curve is less 
[greater] than the angle of inclination @ of the propagation line drawn from thts point. 
Both angles lie between 0 and x and are decreasing functions of s|r]. 


The lemma is obvious in view of (15) and previous results in §1 on the slopes of 
r, s-curves and propagation lines. 


LeMMA 8. If Ty puts a region R of the state plane in (1-1) correspondence with a 
region X of the (x, t)-plane and if the Jacobian J of T., never vanishes in R, the curvature 
of an r|s|-curve in X has a fixed sign and the parts of the propagation lines drawn on the 
convex side do not intersect. 


4. The first initial value problem. Returning to the problem formulated in the in- 
troduction, the correspondents of the initial states (0, po), (0, p2) of the gas are repre- 
sented by the points Po(ro, so), P2(r2, s2) of the state plane in Figure 3a. Both Po, P» 
lie on the line of zero velocity r-+s=0, with rg¢>r2, since po> pr. 

From Lemma 4 we observe that Po lies directly underneath a point Q(ro, s:) of 
the compatibility curve of P: and therefore, according to Lemma 5, a shock and 
buffer wave are generated simultaneously in the gas at x=1. In Figure 3b the shock 
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line from A(1, 0) is AQ” and the buffer region is Pj AQ’. States in the regions 
OAP? , Q'AQ"’, Q’'Ax are represented by points Po, Q, Pe respectively in Figure 3a. 
It is obvious from symmetry considerations that a shock line AQ” and a buffer region 


$ t 























a=0 u=r+s=0 


Fic. 3a. Fic. 3b. 


Pj AQ’ emanate from A(—1, 0) and that states in the regions OAPs , 0’AQO”’, 0’ Ax 
are represented by points Po, 0, P: in Figure 3a with Q the reflection of Q in the line 
r+s=0. In the buffer region emanating from A|[A], we have r=ro[s=so] and the 
equations of the propagation line are 


x — B(ro, s)t = +1, oS SF SS [x —a(r, sot = —1, nsrs ro] 
(ry = —5,). (22) 


As s|r| ranges from so[ro] to si[r:] the propagation line from A[A] turns from 
AP§ |APé | to AQ’[AQ’], with t=Ge" where Go=G(po) at Pé. These propagation 
lines intersect on the t-axis above Py to assign different states to their intersection 
points. We avoid such a physical impossibility by terminating them on the arcs 
Pj Q’, P/Q’ in Figure 3b. The propagation lines assign the states on QPoQ to the 
points of O’PéQ’ and we seek a T, which carries OPQ into Q’Pé Q’ and assigns the 
same states to the same points of the latter arc. A comparison of (7) with (22) leads 
to the following initial value problem. 


THE First INITIAL VALUE PROBLEM. Given two constants ro, So, find a solution w™ 
of (16) for which w(r, so) = —1, w(ro, s) = +1. 

Before giving the solution for the general adiabatic case, we recall a few facts 
concerning the resolvent" of (16). This resolvent is a two parameter family of solu- 
tions v=v(r, 5; fo, So) of the conjugate equation (7 —s)v,.+m(v,—v,) =0 meeting the 
initial conditions »,(7, 50; fo, So) = +1, vs(70, 5; fo, So) = —1, and is given by 


fo—~-s\™ 7F=— % T= fs 
v=(r—- ro ( ) Fact = om; ms — m3 2; : ) 


fe 2a So — To 5S To 


7 — $o\™ S= % S— Se 
- (ss) ( )FaCd = ms m5 — m2 ; ), 


’ 
9g = Se fo — So f= Ze 














M4 See I, in particular §3 and §5. 
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where F; is Appell’s first hypergeometric function of two variables. The solution w\ 
of the first initial value problem is obtained by replacing m by —m and changing the sign 
of the resolvent. 

Monatomic gases present the simplest mathematical problem and from now on 
they will receive our attention exclusively. For them m=1 and 











won Mew +G-me-w A-P--d 4 
r—s r—s 
the last equation holding provided 79+s9 =0. 
Comparison of (17) and (23) yields 
R= ro - r., 5 = s — So (24) 
so that (18’) and (19) become 
r — s)* + 6(rs + ros rs + ros 
en (ry + s) ( ( oe) Gai, i mncthnenens (25) 
(ry — s)* (r — s)? 
{ = O(r — s)-(as — 292) *) = — Or — s)-4*(Br — 27%), (26) 


where the superscripts record that w=w in T,. In the subsonic region of the state 
plane a>0, 8 <0 and, therefore in this region 
€ <4: ¢ >6@ FF >6 (27) 

The square PoQP,0 in Figure 3a is termed the primary region. As po increases 
from p2 the primary region expands from point P? till eventually P,; leaves the state 
plane. We consider only values of po for which the primary region lies entirely in the 
subsonic region and forego examination of the several interesting cases which arise 
when this is not the case. 

Arc Pg Q’, the transform of PoQ by 7,, is tangent to AP¢ at Pj and has slope 
tan 6=1/a>0. Since t increases by (27) and the acute angle @ decreases by Lemma 
7 with increasing s, arc Pj Q’ is concave downwards. Likewise arc Q’P/ the trans- 
form of QP; by 7“), is concave downwards. From (23) and the corollary to Theorem 
1 it follows that arcs Pj Q’, O’P/ are concave downwards. 

The boundary PoQP,0 of the primary region and the boundary P/ Q’P/ 0’ of 
its transform under 7,“ are in one-to-one correspondence with J“ >0 holding in 
the interior of the primary region. It follows!* that the interiors of the two regions 
are in one-to-one correspondence to assign a unique state to each point of the region 
PJ Q’P!0Q’ in Figure 3b. 

5. The second initial value problem. To extend our knowledge of the states of the 
gas we draw propagation lines from the points of the arc Q’P/. These propagation 
lines are tangent to r-curves on Q’P/ and according to Lemma 8 do not intersect on 
the convex side of Q’P/. 

The equation of the propagation line from Q’ is 


rots 


To — Si 





? 


’ (1) 
ot ar ey ee =~ Ps 


8 A beginning in this direction has been made by Rand, loc. cit., 
16 See, for example, G. A. Bliss, Fundamental existence theorems vol. III, Amer. Math. Soc. Col- 


loquium Publications, reprinted 1934, p. 42. 
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and the equation of the shock line from A is, from (13), since re+s,=0, 


(ro + $1) (ro asi $1)* nk 


(ve — 5:)° — (r2 — S2)8 





The two lines intersect in a point Q’’ with coordinates 


1279 1279 = (ro — $1)* — (42 — 52)? 
x”! = 1 4 Pres SPN (ro + $1)(To = $,)’, dd 1 


A(ro, $1) = A(ro, $1) ro — S1 





, (28) 
where 
Ar, s) = (ry — s)*4 — 2(2r + 5)(r2 — Se)’. (29) 


Referring to Lemmas 1 and 2 a solution w™ of (16’) transforms QP, into an s-curve 
propagated from Q’P; and containing Q” if 


we (r, 51) = we (r,s), nSrSry we (ro $1) = 2 — Fro + 2s1)t”, (30) 


the latter condition determining w,‘ uniquely on QP,. 

















$s t 
v=f-S$=0 
r 
\ PR @ 
\ + A=0 
R P, 
Q P, 
=0 u=r+s=0 
Fic. 4a. Fic. 4b. 


Let the arc Q’’P#’ in Figure 4b indicate the prolongation of the shock line AQ’’. 
On the right of Q’’P#’ the state of the gas is P2(re, se) and the states immediately on 
the left of Q’’P/’ are represented by points (r, s) on the compatibility curve (21). 
Thus Q’’P' is the transform by Ty of the compatibility curve QP2. 

On the one hand the slope of Q’’P2’ is 


dt tr’ +t, 3 ter’ + t, 


dx xv + Xx, ae (r + 2s)t.”’ + (2r + s)te 
where r =r(s) is defined implicitly by (21) and its derivative r’ is 
8(r — s)® — 20(r2 — 52)3(r + 2s)? — 3(r2 — 52)§ 


8(r - s)® — 20(r2 — Se)3(2r + 5)? — 3(r2 — S2)§ 
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On the other hand, from (13) 








a = (31) 
dx (r+ s)(r — s)8 
and a comparison of the two results yields the condition 
ts u(r, 5) . : 
—=— r’ = 6, u(r, s) = (r — s)4* + 2(r + 2s)(r2 — 52), (32) 
es X(r, s) 
along QP», or 
i. - i+ eo. +a, =O (32’) 


THE SECOND INITIAL VALUE PROBLEM. To construct a solution w of (16’) meeting 
the conditions (30) on the side QP, of the primary region and the condition (32’) along 
the arc QP» of the compatibility curve. 

From (24) the first condition in (30) is met by taking 


9 


R= r - r, S(s;) = Pi — $0, (33) 


in (17) and, from (28), the second condition determines 
1656 Ke (34) 


The parametric equations of the arc Q’’P/’ are obtained by placing s=s, in (187) 
and substituting for R, S(s:), S’(s:) from (33), (34). In particular it is readily verified 
that x(P}’ )> 1. 


Taking condition (32’) in the form (32), and substituting for ¢t,, ¢, from (19) with 
R=r—r’ it will be found that this condition becomes” 
(7 wa s)? (1) (1) 


. 2+6 1+ 
S” +2——_S' + 6 ——_S = —(t, — st, ), (35) 


r—s (r — s)? 3 


where 6, f{”, ¢? are the rational functions of r, s defined in (32), (26), and r is the 
algebraic function of s defined in (21) with re+s2=0. Thus to obtain the solution w 
of the second initial value problem we set R=r,,—r* in (17) and choose S to be the solution 
of the ordinary differential equation of second order (35), subject to the initial conditions 


in (33), (34). 


LemMA 9. The value of 6 at a point P of the compatibility curve QP2 tends to + 
as P tends to Pz. More precisely 6 is a positive regular analytic function of the parameter 
v on QOP2, except at v=ve, where it has a pole of the third order and a Laurent expansion 


of the form 
§= ve(v = v2) (1 +.--+-), (36) 
To prove 6>0 we have 
A= ov _ 0») _ 3urs, = ov - v) a uve, 


17 The form of the second member in (35) is due to Rand. 





| 





1947] A PROBLEM IN THE PROPAGATION OF SHOCK 341 


and the equation of the compatibility curve QP, 

15v.u'0 = (v - v)(0" - 09). 
It is obvious that p>0 for v>v2, and A\>0 follows from Ay =v?(v* —v3)? —9u28 >0 for 
v>vw.. To establish this inequality, one multiplies the equation of the compatibility 
curve by 3v3 and observes that 3v3(v5 —v3) < 5v5(v? —v}) holds for v>v,. Now r’>0 by 


Lemma 6 and therefore 6>0 for v>w by (32). 
From (20) one has u=r+s=v—v2+2k(v—v2)*+ - +--+, so that 


N= 60x(0— m2) +--+, w= 6n(o— 0) +---, & = (v— 0) (1+ ---), (37) 


hold along QP2, and the Laurent expansion for 6 then follows from (32). 

It is apparent from (20) and Lemma 9 that the coefficients of the differential 
equation (35) present a singular point at s=5». 

Lemma 10. The introduction of v as independent variable in the differential equation 
(35) leads to a differential equation for V =S(s(v)) in which the coefficients are regular 
analytic functions of v for v2 Ve. 


Retaining the prime to denote differentiation with respect to s and indicating 
differentiation with respect to v by a dot, so that S’= V/s, S’’=(sV—Vs)/s*, the 
differential equation (35) becomes 


qa) (1) 


: F - 252 
V¥ + (2(2 + 8)8/v — 3/8 V + 650 (1 +8)V = = gf? af, (38) 


in which the coefficients are regular analytic functions of v for v>v, by Lemmas 6 
and 9. Moreover if the coefficients are expanded in powers of v—v,_ using (20) and 
(36), it will be found that they are also regular about 2». 


LEMMA 11. Provided po—p2>0 is sufficiently small, S and S' are negative for 
Sos Ss, with S tending to a finite limit and S’ to — ~ as s approaches s2. 

Since the coefficients in (38) are regular at v=, the solution determined by 
V(v0) = Vo, V(vo) = Vo may be expanded'® in a power series in v—%2, ¥9 —V2, Vo, Vo pro- 
vided the absolute values of these quantities are sufficiently small. 

Taking v for the value of v corresponding to point Q on the compatibility curve, 
Yo —ve can be made arbitrarily small by taking po— se sufficiently small, with the co- 
ordinates of Q given by 

ro = To + Vo — Ve + K(Vo — v2)? +---, $1 = Se + «(09 — v2)? +---. (39) 


The initial conditions for S in (33), (34) lead to the initial conditions 


2 2 , u(ro, $1) 
V(v) = 51 — fo, V(%) = 25(00)| 5 — 21 we, (40) 
A(ro, $1) 


for V. From (37), (39) we obtain the expansions 


V(v) = — v2(v0 — v2) — (% — m2)? +---, V (v0) = — 2(%— 2) +---, (41) 


18 J. Horn, Gewéhnliche Differentialgleichungen beliebiger Ordnung, Sammlung Schubert, vol. 50+ 
Leipzig, 1905, pp. 27-28. 
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valid for sufficiently small | v9 —v2| . It follows from (38) that the expansion of V(vo) 
in powers of v9 —v2 begins with a term of at least first degree in v9 —v2. 

When the expansions (41) are substituted in the expansion of the solution V in 
powers of v—v2, ¥9—v2, Vo, Vo it appears that V may be expanded in powers of v—12, 
vo—v2 provided |v—v|, |v —v| are sufficiently small. To obtain the linear and quad- 
ratic terms of this expansion, we substitute from (41) in Taylor’s series 

. ie V (v0) 
V = V(v) + V(v0)(v — v) + — (9 — %)* + --- 


to obtain 
V = — 09(00 — 2) — 2(0 — v2)(v — v2) + (00 — m2)? +---, 


the third term in Taylors series being neglected since V(vo) contains the factor v9 —2. 

It follows that both V, V are negative for v,<v Sv for sufficiently small vp —v.>0. 
It is clear that S tends to a finite negative limit as s tends to sz and, since $ is positive 
and tends to zero as v tends to ve, one concludes that SS’ tends to — © as s tends to s», 
provided, of course, that po—/z is sufficiently small. 

The subregion P2QPiRP; in Figure 4a of the primary region is termed the second- 
ary region. 

Lema 12. The partial derivatives t®, t and the Jacobian J = —2G1?t® of Ty 
are negative in the secondary region for sufficiently small po—p2>0. 





We take R, S in (19) as determined by the second initial value problem and find 
(2) 1 (1) 39’ 9S (2) 1 (1) 3 ar 6S’ 9S 
ee Sew Ff ie, a a = ~ oe anne a aeneeney 
2 (r—s)®? (r—s)4 2 2(r—s)? (r—s)® (r —5)4 


from which #2 <0 follows from (27) and Lemma 11 for sufficiently small po—p2>0. 
To prove #! <0 we have 
2+6 1+6 7 iG —7" a oth 


S’ + 6 — S = i, = @. ), 
7—s (# — s)? 3 ( ) 











Ss” +2 


where #=r(s) is the function of s defined in Lemma 6 and 6=4(7, s), £? =1(, s), 
i =1(#, s). When S” is eliminated from ¢@ it is found that t?)=AS’+BS+C, 


where 


6 6 7#—s 
se ane [1+>-——]> ar — yr — 9)*( 


(*# — s)(r — 5s)? 2 r—s 


9 —7—s\* = " 
———_——_—_—_—— [1+3-( <)|> 97 - se - 9+ (5-4), 
(ry — s)?(# — s)? r—sS r; 


1 «) Fe er) ee a ee 
C = — og + 5(- —) (57, 7s "im < - (50, + ig) 


8 
Fs 


2 2 


Nm) ol 
m al 
- lo 
” 





in view of the inequalities 


2ny=n—sy<r—-sSsr—s <9 — So = 2ro, 
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valid in the secondary region. ¢{? has a negative upper bound and ¢" a positive lower 
bound in the secondary region independent of po. Moreover ro, 7; tend to 72 as po 
approaches ps. In view of Lemma 9 we have A>0, B>0, C<0, and therefore 1? <0 
by Lemma 11 for sufficiently small po— pz. 

We shall now investigate the mapping by Ty of the secondary region upon the 
(x, t)-plane. Taking R, S in (18’) as determined above in the solution of the second 


initial value problem, Ty is 








es (r+ s)(3ro _ r) — drs. r+s wy+rs |, 
x 4 = ee —_—— — a = 
(ry — s)? (r — s)? (r — s)? 
ro — 4s 3S 3. Ss 





t@ = 3 — ie amen, eee 
(r—s)*? (r —s5)8 2 (r — s)? 





From Lemma 11 it follows (at least for po— 2 sufficiently small) that x, ¢@ become 
infinite as the point (r, s) of the secondary region approaches the side RP2. We shall 
accordingly consider first the mapping by T.@ of the subregion UQP,TU, the line 
TU being parallel to RP». 

Sides P,Q, TU transform into s-curves P{’Q’’, T’’U’". From Lemma 12 ¢ decreases 
and from (5) x increases as r increases along P:Q, TU. We conclude from Lemma 7 
that P/’Q’’, T’’U” are concave downward as shown in Figure 4b. 

Side P,T transforms into an r-curve P{ T’’ which is concave upwards. 

Arc QU of the compatibility curve transforms into the shock curve Q’’U’’. Along 
QU x and t are monotonic decreasing functions of v, as is the slope dt/dx of Q’’U”, 


for, from (31) 





d {dt Av + ps 
Ace) (¢+s)(r—s)* 
inasmuch as \>0, »>0 hold on QU. The shock curve is accordingly concave upwards 
to imply that the velocity of propagation of the shock decreases as ¢ increases. 
Finally we let U approach P: along QP2. The s-curve T’’U’’ recedes to infinity in 
the (x, #)-plane and the secondary region, exclusive of side RP2, is accordingly mapped 
in (1-1) fashion by JT, upon a region indicated by P?'Q’’P/ R”’ in the (x, ¢)-plane 
to determine the states of the gas in this region. 
The slope of the r-curve P{’ T’’ at T’’ tends to 1/a(r, s2) as T’’ recedes to infinity. 
The slope of the shock curve at U”’ is, from (11), 


dt p2 p— pe 
pate pl 
dx p p— pr 


where p denotes the density for the state U and p=p(p). As U-+P: we have p—:, 
which implies dt/dx—1/G:. This means that the velocity of propagation of the shock 
tends toward the local velocity of sound in the exterior body of gas through which 
the shock travels as it recedes to infinity. 

The determination of the states of the gas in the region Pj’Q’’ P/’ R’’ may now be 
left to symmetry considerations or Theorem 1. 
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6. The third initial value problem. We take up the problem of determining the 
states of the gas in the region of the (x, ¢)-plane lying above the curve R’’P/ P/ P/ R”’ 
in Figure 4b. 














$ t 
v=r-s=0 
r 
LR 
M N P,; A=0 
R 
2 
P, 
a=0 u=r+s=0 
Fic. 5a. 





Propagation lines drawn from P{’R”’ in Figure 5b have slope 1/8 <0 and, from 
Lemma 8 can intersect only on the concave side of P/’ R’’. We shall prove that they 
do not meet in the region x >0 if po—pe is sufficiently small. Since ¢ is a monotonic 
decreasing function of s by Lemma 7, it will be sufficient to prove that the t-intercept 
T of a propagation line is a monotonic decreasing function of s. 

In the equation of a propagation line t=8-'x+T we replace x, t by the coordinates 
of a point on P{’ R”’ obtained from T,,@ to obtain 


us w (rn, s) . (nn, —s)\S’'+S+n— 1 


B(ri, Ss) B(r1, s)(r1 — 5)? 





from which 
dT Br, 5)(11 — 5) S” + 45(71 — 5)S’ + 48S + 4s(ry — 10) 
ds (11, s)(r1 — 5)? 





After S’’ is eliminated by (35) it will be found that dT/ds<0 holds for sufficiently 
small po — p2. The principle of the argument is essentially the same as the one employed 
to prove that £” <0 in Lemma 12 and is omitted. 

From symmetry considerations propagation lines drawn from P/’ R’’ do not in- 
tersect in the region x <0. Propagation lines drawn from P{’ R"’ and P{’ R”’ symmetri- 
cally placed with respect to the t-axis intersect upon it and, excepting the two drawn 
from P/’, P/’, assign different states to their points of intersection. This is avoided 
in Figure 5b by terminating the propagation lines on arcs P{’’N’’’, P{” N’’’, the 
coordinates of P/’’ being x =0, t= —w®”)(r,, s:)/B(r1, 51). 


EEE 
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By Lemma 1 an r-curve, the transform of P:R by Ty will be propagated from 
P{' R”’ provided 
(3) (2) 
Ww, (71,5) = We (1,5) for se S55 51, (42) 
and by Lemma 2 will contain P/ ’, in case 


(3) a(ri, $1) @) (2) 
w, (fr, $1) = — WwW, (71, $1) = — Ws (F1, 5). (43) 


B(r1, $1) 

At points symmetric to the t-axis states have the same density and opposite veloci- 
ties. From the corollary to Theorem 1 this will be the case in the region above 
N’"'P! 'N’”’ if this region is the transform by 7, of a region in the state plane sym- 
metric to the line r+s=0, provided w(r, s)=w(—s, —r) holds in this region. 





THE THIRD INITIAL VALUE PROBLEM. To construct a solution w® of (16’) meeting 
the symmetry condition w®(r, s)=w®(—s, —r) and the initial conditions (42), (43) 
on the side P; R of the secondary region. 


The solution of this initial value problem 


S(— r) + S(s) 
wr, 5) = — ——__ (44) 
r—s 
is obtained by setting R= —S(—r) in (17), where S(s) is the function entering in the 
solution of the second initial value problem. 

The symmetry condition is obviously fulfilled. From (33) we find w®(n, s) 
=w(r,, s) and (42) follows by differentiation. Condition (43) is likewise a conse- 
quence of (33). 

The subregion P:RP2RP, of the secondary region in Figure 5a is termed the 
tertiary region. 

The mapping by 7. of the tertiary region upon the (x, ¢)-plane is not (1-1). 
If R is replaced by —S(—r) in (19) one obtains 

3 (r — s)2S’"(— r) + Ar — s)[28’(— r) + S"(s)] + 6[S(— r) + S(s)] (45) 
» (r — s)4 





In particular on r= —s (along P,P») 
3 
{@) = geals*S’(5) — 3sS'(s) + 3S(s)]. 
, s ; 


For po—p2 sufficiently small ¢ is positive along P,P: in view of (27), (35) and Lem- 
mas 9 and 11. On the other: hand, if we fix’r in (45) and allow s to approach 5¢ it ap- 
pears that /* eventually becomes negative because of the behavior of S, S’ as s tends 
to se. Hence there exists a subregion P;P:MP, of the tertiary region within which 
t? is nositive, except along MP2 where { =0. 

Differentiating equations (5) partially with respect to r and s and eliminating x;., 
we find that ¢ satisfies the partial differential equation (a—8)t,, =8.t,—a,t, which re- 
duces, in the adiabatic case, to 
Se nm. 
(7 — st, = — —— 


-(t, — 6). 
, oa" ) 
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PLATE 1, Variation of density p with distance x for fixed time ¢. 
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PLATE 2. Variation of density p with time ¢ for fixed distance x. 
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It follows that in a region in which f, is positive and ¢, is negative, /,(¢,) is a monotonic 
increasing function of s(r) for fixed r(s). 

From t®(r, s) =t®(—s, —r) we have ¢®(r, s) = —#(—s, —r) and thus 4 = —1° 
along P,P». It follows that ¢ is negative along P,P2, and consequently is negative in 
a suitably restricted region containing P:P2 When we procede from P,P» to the left 
along a line s=const., ¢ decreases and accordingly is negative everywhere in 
P,P.MP,,. 

Along MP; dt\*/ds is positive and J = —2Gt changes sign as (r, s) crosses 
MP». Thus the mapping of the tertiary region cannot be (1-1). 

It is clear that ¢ vanishes along an arc MPs, the reflection of MP2 in P,P». 
Within the region P; MP.MP, we have J® >0. The application of T,,@ to the tertiary 
region to obtain further information about the states of the gas must be restricted to 
a square P,NP;N within the region P,MP.MP, with sides on P:M, PiM. Such a 
square is carried by T,@ into the region P{/’N’’’P{’"N’"’ of Figure 5b, within which 
the states of the gas may be regarded as known. 

The prolongation of the solution into the rest of the (x, ¢)-plane still awaits solu- 
tion and it should be noted here that further extension of the solution may modify 
the states assigned above to the region R’’P{'Q’’P?’ and its reflection in the ¢-axis. 

7. Graphical presentation of variation of density. Plates 1 and 2 portray the varia- 
tion of density (and thus the pressure) of a monatomic gas for sufficiently small 
Po— 2 in so far as our analysis permits. They were obtained by comparing Figures 5a 
and 5b with the aid of (2). No attempt was made to indicate quantitative changes in 
the density, or to determine the curvature of the curved portions of the graphs. The 
small circles indicate points at which pz and p,; undergo jumps. Figures VIII, LX of 
Plate 2 are based on the conjecture that x(Q’’) <x(N”’). 

The coordinates of various points in Figure 5b have been computed from formulas 
found in the text and are given below. Once po, pz are given, fo, 72 are determined from 
(2) and s, may be determined graphically as the s-coordinate of Q in Figure Sa. 
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THE PROPAGATION OF A SPHERICAL OR A CYLINDRICAL 
WAVE OF FINITE AMPLITUDE AND THE PRODUCTION 
OF SHOCK WAVES* 


BY 


YUNG-HUAI KUO 
California Institute of Technology 


1. Introduction. When a mass of gas is set into motion by a sudden rise of pres- 
sure which possesses either a cylindrical symmetry or a spherical symmetry in the 
case of an explosion, pressure or density will be propagated into space as a cylindrical 
or spherical wave of finite amplitude in a manner different from that of the propaga- 
tion of sound. The most conspicuous phenomenon of such a non-linear wave motion 
is perhaps the appearance of a shock wave. In the case of plane waves of finite ampli- 
tude, the problem was studied independently by B. Riemann! and S. Earnshaw.? It 
was shown that when a compressed slab of gas is released, two progressive waves are 
produced travelling in opposite directions, with constant deformation in the wave- 
form during the course of the propagation. Eventually both waves develop into shock 
waves. 

With regard to the spherical or cylindrical compression waves, the situation is 
quite different because the amplitude of the wave falls off at a much greater rate than 
for plane waves, while the wave propagates from the center of disturbances. The 
question is whether this rapid diminution of amplitude would prevent the formation 
of a shock. J. J. Unwin’ has calculated a specific example of motion produced by a 
sudden release of a compressed sphere of air, and concluded that there is no indication 
of the development of a shock wave. Inasmuch as he adopted a numerical method for 
one special case, the concludion reached cannot be regarded as general. In fact, 
W. Hantzsche and H. Wendt‘ considered a similar problem, where the sphere had a 
finite radius and expanded with the speed of sound into still air. The motion, in its 
early stage, is supposed to be continuous in pressure or density and velocity. But after 
a finite duration, the wave-front becomes a discontinuity surface characterized by an 
infinite velocity gradient in spite of the diminution of amplitude. 

In view of these disagreeing results, it is felt that it is desirable to investigate this 
problem from a broad standpoint taking account of all initial boundary conditions. 
The problem of explosion such as the burst of a bomb is only one of many similar 
problems and, to be sure, the most interesting one. According to G. I. Taylor, the 
physical process taking place during an explosion can be treated, as a combination of 
two problems. The first problem is concerned with the effects produced in the atmos- 





* Received May 15, 1946. 
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phere by a rapidly expanding spherical or cylindrical solid shell which compresses the 
surrounding air. In this case the motion of air in contact with the shell is completely 
prescribed by the motion of shell itself. The second problem deals with the motion 
produced by a compressed sphere or cylinder of air which is suddenly released. Each 
one of these constitutes a separate mathematical problem. To enlarge the scope of 
this discussion, the very meaning of the term explosion will be understood here as 
any process that is capable to create a pressure disturbance with spherical or cylindri- 
cal symmetry, propagating as a wave of finite amplitude. 

An explosion is assumed to take place, during a short interval of time, in an in- 
finite space which is filled only with air not abstructed by any solid bodies. Since the 
coefficients of viscosity and heat conduction for gases are generally very small, so 
long as the motion is continuous, the air may be regarded as non-viscous and non- 
conducting. The thermodynamic change of state of a fluid-particle along the path is 
then adiabatic; and if, initially, the entropy of the air is uniform throughout the space, 
the motion is isentropic. For the first problem stated above this condition is satisfied. 
Namely, at the moment the shell starts to expand, the outside air may certainly be 
assumed to be at the standard conditions. After the shell has started to expand, it 
compresses the air and, thereby, sets it into motion; but, during this process, no heat 
has been imparted to the air, its thermodynamic state must remain on the same 
adiabatic curve. In the case of a compressed sphere or cylinder of air, it is reasonable 
to assume that the pressure or density was built under adiabatic compression at all 
points. Hence as long as the motion is continuous, it will be isentropic. 

The present study reveals that such a continuous and isentropic motion generally 
does not exist in the whole field. This type of motion breaks down when a “limiting 
line” appears, which would make the solution multi-valued. This would be impossible 
unless the motion is discontinuous. Hence, the appearance of a “limiting line” serves 
to indicate the necessity of presence of a shock wave in the actual motion. After the 
shock is formed, the Rankine-Hugoniot theory asserts that the process through which 
a fluid-particle has undergone by crossing the shock-front is irreversible and, conse- 
quently, the entropy increases in a discontinuous manner. The jump in entropy is 
not constant, however. It varies as the shock wave propagates, because the conditions 
at the shock change with time. As a result the motion behind such a non-uniform 
shock cannot be isentropic. Therefore once the “limiting line” appears, isentropic flow 
cannot be maintained and the resultant flow cannot be analyzed by the present 
method. 

The mathematical condition for the appearance of a “limiting line” in the case 
of a spherical or cylindrical isentropic motion is that one of the two families of char- 
acteristics admits an envelope, just as in the case of a plane wave. Along this envelope 
the accelerations of the fluid-particles are infinite. In fact, a closer examination in- 
dicates that the motion generally must break down even before the “limiting line” 
is reached. It then seems that any motion of a compressible fluid has a tendency to 
develop a shock wave and that the effect of the “spreading” in the case of a non-linear 
spherical or cylindrical wave plays but a minor role. 

2. Differential equations of motion. The motion under consideration is supposed 
to be axially or spherically symmetric, i.e., at any instant the velocity u, pressure p 
and density p depend on the time and the radial distance x only. If the effects of vis- 
cosity and of body force are neglected, the equations governing the motion are 


—EE 
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Pe are er (2.1) 
p 
au 
pet Upz + Au. + =) = 0. (2.2) 
x 


Here the subscripts denote the partial derivatives with respect to the variable indi- 
cated by the subscript; a=1 for a cylindrical and a=2 for a spherical wave. In each 
case, the variable x will be interpreted differently. Furthermore, it is assumed that the 
motion is continuous and that the effects of viscosity and heat-transfer in the fluid 
can be ignored. If initially constant, throughout the fluid, the entropy then remains 
constant. In other words, for an ideal gas the relation between the pressure and den- 
sity is 


p = Kp’, (2.3) 


where y stands for the ratio of the specific heats and K is a constant. With a set of 
appropriate initial conditions the mathematical problem can then be solved, at least 
theoretically. However, we may understand the singular behavior of such a solu- 
tion and the conditions for its existence without actually solving the differential 
equations. 

By eliminating the pressure with the aid of Eq. (2.3) and by introducing the square 
of the sonic speed as a variable in the place of the density, we reduce Eqs. (2.1) and 
(2.2) to 


u,+ uu, +, = 0, (2.4) 


au 
VU: + Uv: + ow. + wn) 
x 


0, (2.5) 


where 
fo=c, B=ay-l, 


and c is the speed of sound defined by Vy(p/p). This system of differential equations 
is of the hyperbolic type, the two families of real characteristics C being determined by 


(dx — udt)? — Bodi*? = 0, (2.6) 


where v is positive. 

As it stands, this system of equations can reveal but little information concerning 
the behavior of the solution. To expose such properties, one has to transform the dif- 
ferential equations to a new coordinate-system and then study the condition under 
which the transformation would be valid. In the case of a steady irrotational motion, 
this is well-known as the hodograph method which has been effectively and success- 
fully applied by W. Tollmien® and H. S. Tsien® in investigating the two dimensional 
and three dimensional isentropic motion respectively. By a slight modification, it 
can also be applied to the present problem. To this end, the following one-one point- 
transformation is introduced 





5 Tollmien, W., Grenzlinien adiabatischer Potentialstroémungen, Z. angew. Math. Mech. 21, 140 (1941). 
6 Tsien, H.S., The “limiting line” in mixed subsonic and supersonic flows of compressible fluids, N.A.C.A. 
Tech. Note 961 (1945). 








352 YUNG-HUAI KUO [Vol. IV, No. 4 


u = u(t, x), v = v(t, x). (2.7) 
We have 
y by 
“=-») te = =e 
J J 
a bes 
Vr —_ = ae v7, =-——»? 
J J 


provided the Jacobian J(u, v) =t,x, —t,x,#0. Equations (2.4) and (2.5) will then be 
transformed into 
0, (2.8) 


xy — uty + tu 


0. (2.9) 


aBuv 
Ly — uty + Buty — —— (but, — bye) 
x 


This system of equations can be simplified considerably by introducing a function 
x(u, v) defined by 


x— ut = Xx; i= — Xp, (2.10) 
so that Eq. (2.8) is satisfied identically while Eq. (2.9) reduces to 


aBuv 2 
Xuu ~~ BUX vv a Ca ee OP ie XX vv) = Xs 
x 


(2.11) 


The corresponding characteristics [' in the u, v-plane are determined by 


aBuv 2aBuv 1 Xuu — Xv 
(1 — xs) dy? — Xurdudv — aBuv (— + ) aw = 0. (2.12) 
x 


x au x 











3. Limiting line. The relationship between the characteristics C and I associated 
respectively with the differential equation in the ¢, x- and u, v-planes has an important 
bearing on the singular character of the solution and its elucidation often contributes 
much toward the understanding of the nature of the physical problem. For this pur- 
pose, we first transform the differential equation (2.6) by means of the following pair 


of relations: 

dx = (Xuu — Xv — UXuv)du + (xuv — UXvr) dd, 

dt = — Xurdu — xrrd2. 
Substituting in Eq. (2.6) together with Eq. (2.11), we bring the equation of the char- 
acteristics C into the form 


Buy 2aBuv 1 Xuu — Xv 
J (1 uneocees x) dv? — —— yxy,dudv — ad (— + =) int =Q. (3.1) 


x x au x 


This shows that if J#0, the characteristics C in the ¢, x-plane correspond to the 
characteristics [ in the u, v-plane. However, circumstances may arise such that 


2 
J(u, v) = XuuXeov ~ Xuv ~ XvXvv = 0, (3.2) 
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while 


aBuv 2aBuv 1 Xuu — v 
(1 —- x) dv? — ——— xxu,dudv — aBuv (— + x2) du? ~ 0 


x x au x 





and the characteristic equation (2.6) is again satisfied. This means that if a point 
moves along a line \ defined by Eq. (3.2), the corresponding point will describe a 
line / in the ¢, x-plane, having the same tantents as the characteristics C. It does not 
coincide, however, with any one of the characteristics C. This may be proved as 


follows. 
The differential equation for the path s of a fluid-particle in the ¢, x-plane is 


(= os 
=) = 4. (3.3) 


The corresponding path o in the u, v-plane is given by 


dv Xuu — Xv 
(=) -=-— ". (3.4) 
du}, — 


Now the differential equation for one family of characteristics, say T’,, is 


dv Ree “ Xe + / BY Xuv 

(2) = ent Bae ids 
du ry Xuv + / By Xvv 

On the other hand, the vanishing of the Jacobian, when combined with Eq. (2.11), 


can be written as 








Chie = \/ BY Xvv) (Xuv + / BP Xr) = 0. (3.6) 


(2) (Gm : 
du - xh » 1 Xuv = BY Xvv (3. ) 


(=) ( =) ‘e 
— 2s, ’ 1 uv 
du} , du/r_ 7 


The condition under which this result holds is both necessary and sufficient. This 
shows that the lines \, and \_ are respectively the locus of the points of tangency of 
the path o with I’, and o with T_. Furthermore, the paths ¢ do not have an envelope 
and that of T is 


It is easy to see that 


or 


as VB Xov- (3.8) 


po = 0 


which corresponds to p=0 and is, of course, uninteresting. Hence, it cannot belong 
to either family of the characteristics [. The only alternative is that it is an envelope 
of one family of the characteristics C in the ¢, x-plane. By analogy with the steady 
irrotational motion it is again called “limiting line,” the justification will be found in 


the following section. 
4. The properties of the “limiting line.” Being the envelope of one family of real 
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characteristics in the ¢t, x-plane, the “limiting line” will be entirely in the field of mo- 
tion. It is, therefore, paramount to investigate the behavior of the solution along 


this line. 
Consider first the line element of a path s of a fluid-particle at the “limiting line” /. 


Generally, for any line element one obtains from Eq. (2.10) 
dx = (Xuu — UXuv — Xv)dUu + (xuv — UXvv)d2, 
dt = — xu du — xvrd0. 
Along a path s given by dx/dt=u, we have 
(xuun — Xv)du + xu dv = 0. 
Using this relation to eliminate dv from dx and dt and by regarding u as a parameter, 


we obtain the following parametric equations for the path s: 


XuuXvv ~ 7. — XvXvv 











dx=u dn, (4.1) 
Xur 
XuuXev ~ x. ee XvXvv 
dt = ——_—_——— dy. (4.2) 
Xuv 


According to our previous findings, J =0 yields two lines \, and X_, each of which 
associates with only one group of characteristics T' in the u, v-plane. This shows that 
on the “limiting line” dx and dt both become differentials of higher order and will 
change sign on crossing the line \. This agrees, of course, with the cuspidal nature of 
the singularity. 

Dividing both sides by dx and dt respectively, we obtain the following expressions 


for the derivatives u, and u, along s: 





Xur 
sia EN ONS Ty (4.3) 
UXuuXre = pa XeXee) 
. Xuv 
(t+)s a en ———————— : (4.4) 
XuuXvv ~~ ee — XoXve 


Thus on the “limiting line” the acceleration of a fluid-particle becomes infinite as xu» 
is finite there. This implies also an infinite pressure gradient [see Eq. (2.1) ]. 

The physical state to which J(u, v) =0 corresponds can be readily deduced. It 
can be summarized in the statement that if the Jacobian vanishes, then the motion 
in the immediate neighborhood of the line J=0 is a compressive one. To prove this, 
let us consider the ratios v,/u2, U./vz, u,/u, and v,/v, which, according to the relations 
obtined in Section 2, equal 











UV Xuv afuv J Ut ow 

—> = — Bu-— 4 <> qmmma'¢ ammee 9 —_— = —-u+—) 

Uz Xvv x Xovv Uz Xeov 

Ut Xvv U% Xv» aBuyv J 
—- ey —— — I, — = —u+ po + — 
Vz Xuv Vz Xuv x Xuv 


In the u, v-plane, the expressions on the right-hand side are everywhere continuous. 
At the line \, corresponding to Xu» =VBuxe», they become 
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Ve u Ut u 
—=-c{[1——) <0, —=c{1-——]>0, 
Ur c Uz c 
Uz u Ve u 
ae See —=#cl1—-—)>0. 
Vz c Vz c 


By continuity, the relative signs of the differential quotients hold in the neighborhood 
of the “limiting line.” Thus, we conclude that either v,>0, v.>0 and u,<0, u,<0 
or v, <0, v, <0 and u,>0, u,>0. The first case is exactly the condition for a compres- 
sive motion. Whereas the second case may either correspond to a rarefaction or to a 
change of sign of the Jacobian J(u, v). As the rarefaction does not conform to the 
geometric properties of J =0, the second case corresponds to the second branch of the 
solution and hence can be disregarded. 

5. Lost solution. In the previous sections, we assume that the Jacobian J(u, v) 
does not vanish. Thus the one-to-one correspondence between the ¢, x- and u, v-planes 
is assured and the condition J =0 is restricted to the singular line /. In a special case 
the Jacobian may vanish identically, however. This vanishing of the Jacobian estab- 
lishes a relation between v and uw in the u, v-plane and, as a result, yields a class of 
solution not contained in the transformation (2.7). To study this form of solution, 
let us first set 





v = v(). (5.1) 
The differential equations (2.4) and (2.5) can then be rewritten as 
dv 
uy + (u + -) u, = 0, (5.2) 
du 
dv dv aBuv 
w+ (u ~ + 6) t= — (5.3) 
du du x 


This type of solution has been discussed by K. Bechert’ whose main result was as 
follows. By eliminating x and ¢ the system of Eqs. (5.2) and (5.3) can be reduced to a 
second order non-linear total differential equation, based on the existence of a linear 
relation between ¢ and x. By a slightly different procedure it can be shown that in- 
stead of a second order differential equation one can obtain a first order one of Abel’s 
type being amenable to numerical integration. The main feature of the solution, how- 
ever, can be discussed in the following manner. 
Along u=const., i.e., along 


du = u,dx + u,dt = 0, 


the slope of the curve “=const. equals 


dx Ut dv 
( ) = - ante, (5.4) 


dt u Us du 


on account of Eq. (5.2). Since dv/du is a function of u alone, on u=const. (dv/du), is 


7 Bechert, K., Uber die Ausbreitung von Zylinder- und Kugelwellen in retbungsfreien Gasen und 
Flussigkeiten, Ann. Phys. (5) 39, 169 (1941). 
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constant. Therefore, the curve u=const. is a straight line in the ¢, x-plane. In con- 

formity to the assumption (5.1), there exists a parameter & defined by 
x 

co(t + to) 


where ¢o is the speed of sound at «=0, and ¢ a suitable constant. It is clear that 
£=const. corresponds to u=const. In other words, both v and u may be regarded as 


t= (S.5) 


functions of &. 
If the determinant v’*—6v+0, uz and u; can be expressed in terms of u. We have 








aBuv 1 
“4, = Semone (5.6) 
x vw? — Bo 
apuyv u+ov' 
“= — ———— (5.7) 


x v2? — By 


where the prime denotes the total differentiation with respect to u. Like in the gen- 
eral case, here again the solution possesses a singular line on which the partial deriva- 
tives generally become infinite. Its other properties will be studied presently. From 


Eq. (5.4) it is found that 
dx 
(2)+0 
dt/., 


dt/¢ 
On the other hand, where the singular line X, i.e. the line 


v’? — Bv = 0, (5.8) 


while the characteristics are 


u + »/Bo. 


intersects the integral-curve v(u), we have 


(=) i (=) (5.9) 
——e moe oe 8 5.9) 
dt/. ¥ dt/¢ 


This shows that at the singular point of the solution v(u), the «=const. line be- 
comes the envelope of one family of characteristics C. Hence the envelope is a straight 
line. Furthermore, according to Eqs. (4.1) and (4.2) the parametric equations of the 


path s are 


dx = — ——(v' + \/Bo)(0’ — V/Ba)du, (5.10) 
aBvv’ 

dt = — ———(v + o/fn)(0’ — V/Ba)du. (5.11) 
aBuvo’ 


Since each factor on the right-hand side corresponds to a group of the characteristics 
C, on crossing the line A, where this factor vanishes, the elements dx and dt change 
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their signs. This proves that the line /, the image of \, possesses all the characteristics 
of a “limiting line.” 

It is interesting to note the difference between plane and spherical waves. In the 
former case, Eq. (5.8) would be satisfied identically. This lets the lines u=const. de- 
generate into the characteristics. Indeed, it is also possible for one family of the char- 
acteristics which are straight lines to have an envelope; the differential quotients uz, u; 
are finite, however. Consequently, we have no “limiting line,” in the strict sense. 
This does not mean, of course, that the solution is regular. As a matter of fact, the 
solution already becomes many-valued before this line is reached. 

6. Lost solution: a special problem. From the foregoing conclusions, a compres- 
sive spherical or cylindrical wave always becomes indeterminate when a singular line 
is reached. As an illustration the following special problem is considered. 

Suppose there is a divergent spherical or cylindrical wave propagating with veloc- 
ity Co into still air. On the wave-front, where the motion agrees with the outside con- 
ditions, the state-variables p, p become equal to those of the still air and the velocity 
is zero. The path of the wave-front is then described by 


x =.Co(t + bo). (6.1) 
The mathematical problem can thus be formulated in the following way: 
u=0, when x 2 ¢Co(t + bt), 
ux0O, when x < oo(t+ at 


(6.2) 


A particularly simple case will be the one where both the pressure and the velocity 
are propagated with constant speed. In other words, these quantities depend only 


on a common parameter. 
To simplify the amount of mathematical work involved, the differential equations 


(2.4) and (2.5) will be put into the following equivalent form: 


2c*hz 


(2 — $2)b22 — 26 — Gu + —— = 0 (6.3) 
x 
by introducing a potential-function @(¢, x): 
9 2 B 2 
u = dz, et oee= — Boi. (6.4) 


In the case of a lost solution, there exists a parameter £ defined by (5.5) such that 
£=1 corresponds to the initial curve (6.1). Then, 


b(t, x) = co(t + to) f(é) (6.5) 

and hence 
u(t, x) = cof’(é), (6.6) 
= ae - 5 yn — Bf - |, (6.7) 


where the prime indicates the total differentiation with respect to £, and the function 
f(&) satisfies 
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[c? — co(f’ — &)?]Ef” + 2cif’ = 0 (6.8) 
subject to the initial conditions 
AH =0, fi) =o. (6.9) 


The first condition, namely f(1) =0, is necessary to make c=co on £=1. When the 
conditions (6.9) are substituted in Eq. (6.8), it appears that f’’(1) is arbitrary. We 
need not be alarmed by this situation, but recall that in this particular type of initial 
value problem, the “support” is a characteristic. Physically, this means that the ini- 
tial conditions prescribed in this manner do not “know” the internal structure of the 
motion, because they propagate ahead with larger speed. It is only natural, then, that 
such an arbitrariness should arise which enables us to fit properly the physical condi- 
tions specified. This arbitrariness is only a partial one, however, since for a compres- 
sive motion the sign of f’’(1) is necessarily negative; for on —=1 


(p1)1 + po(uz): = 0, 


according to Eq. (2.2). In a compressive motion (p;)1>0, it follows that 
co, 
(wz). = —f"(1) < 0. (6.10) 
x 


Thus, for any compressive motion the absolute value of f’’(1) is determined in con- 
sistence with the physical process. 

The differential equation (6.8) which determines the interior motion of a mass of 
air, has two singular points in the &, f-plane given by the vanishing of the coefficient 
of f’(£). The geometrical interpretation is evident, when (6.8) is written as 


(c + u — coé)(c — u + Cok) 
dx dx dx dx 
=-|(—) -(— —) -(— (6.11) 
dt C+ dt E dt 0.2 dt g 
that is, when one family of characteristics become tangent to a line £=const., an 
infinite curvature would occur if u is finite there. According to what has been said 
in the last section, this characterizes the “limiting line” of the solution. 


Let us push the discussion a step further. For this purpose only the first order 
terms need be retained. Taking 8 as a small parameter, one has accordingly 


HE) = fol + BAl— +---. (6.12) 
Substituting in Eq. (6.8) we obtain 
[1 — (fo — £)*]Efo' + 2fe = 0. (6.13) 
This equation is free from fo; letting w=f’ we find 
dw w 2 


de & (w—8?-1 


IA 
wr 
IlA 


(6.14) 


Aside from the two singular lines 


w=§+1, (6.15) w=t—1, (6. 16) 
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where the slope of w is infinite, there are two additional singularities (1, 0) and (0, 0) 
where the slope is indeterminate. The point (1, 0) acts as a sort of nodal point which 
makes the initial condition insufficient. The point (0, 0) is a saddle point as locally 
the equation behaves like 

dw 2w 
-—-») (6.17) 
dé E 


which form is obtained by neglecting (w—&)* as compared with 1. 

The situation can now be summarized. The integral curve starting from (1, 0) 
rises as — decreases and eventually intersects with the line (6.15) where it will have a 
vertical tangent at §<1. After it crosses this line its slope changes sign. This causes 
the curve to bend backward again. Thus, & is seen to assume a minimum value. Owing 
to the fact that the origin is a saddle point, no integral curve could possibly cross 
the line £=0. This fact makes the continuation of the solution as far as £=0 impossi- 
ble. 

7. Continuation of the solution. The results obtained in the previous sections 
show that, in the case of the propagation of a spherical or cylindrical wave, a continu- 
ous solution does not exist throughout the domain considered and can be constructed, 
at most, as far as a singular line / in the ¢t, x-plane from a suitably chosen initial data. 
The line / thus acts as a sort of “frontier” into which no solution can enter and at 
which the solution is turned back as a second branch. The domain then is doubly 
covered. Physically, this is impossible and hence must be rejected as a solution. The 
question is: is it possible to connect it with a different solution beyond this line? 

First, consider the line \ as a “support” with a given set of initial data and then 
solve the initial value problem® for a Monge-Ampére equation. Regarding \ as a 
parameter, we have along the line A 





d du dv 
“ = Xe he ties ae (7.1) 
d du dv 
£0 gg # (7.2) 
and hence 
>. au d d 
(XuuXve _ Xu) a = Xvv D i aa Xv- 


Substituting this into Eq. (2.11) we obtain a linear relation between the partial de- 
rivatives: 


aBuv (dx,,/dd | ] aBuv dx,/dr 
uu + a a ae v —_— i —. oo eee: ee ae (7.3) 
, | x { du/ dd " ad x  du/dr 


Since \ is not a characteristic, Eqs. (7.1), (7.2) and (7.3) are sufficient for a unique 
determination of xuu, Xu» and x»»; and consequently a unique integral surface. The 
uniqueness of the solution is sufficient to show that the solution, when transformed 





8 Courant, R. and Hilbert, D., Methoden der math. Physik, vol. 2, J. Springer, Berlin, 1937, p. 344. 
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back to the #, x-plane, will correspond to the very one that doubles back at the “limit- 
ing line.” A continuous solution is thus out of the question. 

The alternative procedure would be to continue it by joining it smoothly at the 
line \ to the lost solution. This is also impossible. Indeed, if this were possible, the 
line \ would have to coincide with the integral curve v(u) in order to provide a con- 
tinuous solution. This is contradictory, because it is easy to show that the line \ does 
not satisfy the differential equation for v(x). 

The other possibility which remains to be investigated is to identify the “limiting 
line” as a shock wave so as to construct a discontinuous solution. This would require 
the continued solution to satisfy the shock conditions. Since, in general, the “limiting 
line” ] is curved, as a result there would be a non-uniform shock wave in the motion, 
for which both the speed and the strength are no longer constant and therefore the 
entropy would be constantly changing across the shock. This very fact makes the 
original assumption untenable. Hence to continue discontinuously a solution with 
entropy constant everywhere is also impossible. 

The problem might be solved, however, if the original hypothesis of isentropic 
motion is abandoned. To include the possibility that a shock wave may exist within 
the motion, the continued solution must satisfy the following more general set of 


equations: 


ae ee. == 0, (7.4) 
p 
au 
n+ ues + o(u.+ —) = 0, (7.5) 
x 
(pe~”)e + u(pp-7), = 0. (7.6) 


The task then is to construct a solution which should satisfy both the initial and the 
shock conditions in a region bounded by the initial curve, the shock line and a char- 
acteristic drawn to the initial curve through the point where the envelope first ap- 
pears. The shock line, however, is not given, it should be chosen in such a way that it 
yields a solution fulfilling all the prescribed conditions. The mathematical problem 
thus turns out to be extremely difficult. 

The author wishes to thank Dr. H. S. Tsien for his invaluable discussions and 


criticism. 
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ON PROJECTILES OF MINIMUM WAVE DRAG* 


BY 


WILLIAM R. SEARS 
Cornell University** 


1. Introduction. The wave resistance of slender bodies of revolution in symmetri- 
cal supersonic flow was calculated approximately by von Karm4n,! by means of a 
distribution of singularities along the axis of the projectile. The individual singularity 
is characterized by a potential of the form ¢:(x, r) = { (x —£,)?—a’r?}—1/2, where x, r 
are cylindrical coordinates, x being measured downstream from the nose of the projec- 
tile and r radially from the axis, £; is the value of x corresponding to the singularity, 
a is the cotangent of the Mach angle of the undisturbed flow, so that 


a = V(U/a)? — 1, 


U and a being the stream velocity and the velocity of sound in the undisturbed flow. 
It will readily be verified that ¢;(x, r) is a solution of the linearized potential equation 
for supersonic flow with axial symmetry 
U? a 3% 1 a 
= +——. (1) 


ee. ax? Or? r Or 


Von Karman calculated the wave resistance by integrating the transport of mo- 


mentum across a cylindrical surface enclosing the body. In his approximation, the 
integral is independent of r and can be evaluated in the limit r-0. The result ist 


R=-af f ss'@ tog | x—t| dvds, (2) 


where R is the wave resistance and f(x) is the function specifying the distribution of 

singularities along the x axis. For bodies of finite length /, f(x) is found to be indenti- 

cally zero for x >1; hence both integrals in (2) can be replaced by integrals from 0 to /. 
For slender bodies, von Karman showed that approximately 


fiz) =~ => (3) 


where S is the cross-sectional area of the body. 

In the present paper we shall amplify the analogy, already mentioned by von 
K4rmA4n, between the wave resistance of a slender projectile and the induced drag of a 
wing. It will be shown that this analogy suggests a useful form for the calculation of 


* Received June 11, 1946. 

** This work was undertaken while the author was employed by Northrop Aircraft, Inc. 

1 Th. de Kérman, The problem of resistance in compressible fluids, Atti del V Convegno della “Fonda- 
zione Alessandro Volta,” Rome, 1935, pp. 222-276. 

+ Von K4rm4n!, Eq. (9.12). It might be mentioned that this formula is most easily obtained from 
Eq. (9.11) of the same reference by first integrating by parts with respect to x in order to obtain a form 
symmetrical in x and £; it will then be found that a double integral carried over half the first quadrant 
of an x, £ plane can be identified with half the same integral carried over the entire quadrant. 
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the wave drag. The properties of projectiles of minimum wave drag for given length 
and volume, and for given length and caliber, will then be investigated. 

2. The induced-drag analogy. Formula (2) for the wave drag can be written in 
the form 


R=- of f'(*)F(x)dx (2’) 


or, after integration by parts, assuming f(0) =f(/) =0; i.e. that the body has sharp 
points at front and rear, 


l 
R= pf s(2)P"(a)dx, (2”” 
0 
where 


F(x) = f f'(&) log | x— £ | dé. (4) 
0 


In the form (2’’), von Karmén’s analogy between the wave resistance and the 
induced drag of a finite wing in the Prandtl lifting-line theory? is evident: f(x) is 
proportional to the circulation distribution over the span of the wing, F’(x) is the 
corresponding downwash distribution, and R is the induced drag. 

It is also useful to put (4) in another form, sometimes more convenient for calcula- 
tion. Let us introduce the coordinates 6 and # defined by 


ll 


l 
=< + cos @), 058, 


; (5) 
cm 7 tt + oe, ver 
The expression for F(x) then becomes 
l ° 
F(x) = —f t’(&) log | cos @ — cos 3 | sin ddd, (6) 
0 


provided that /{f()d~=0, as is always the case for closed bodies, in accordance with 
(3). Now the definitions in (5) can be taken to cover the range —7 S0S7, —7 Sd Sr, 
and f’(£) can arbitrarily be defined to be an odd function of 3. Then (6) can easily be 


put into the form 
9-9? 
2 


sin ddd 








l F 
F(x) = —{ f'(&) log sin 


or, after integration by parts, 


6-38 
dé. (7) 





1 . 
F(x) = — =f cot 





2L. Prandtl, Tragfligeltheorie I, from Vier Abhandlungen zur Hydrodynamik and Aerodynamtk, 
Gottingen, 1927, pp. 9-35. 
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The induced-drag analogy pointed out above suggests that f(£) be expanded in a 
sine series; this is the usual technique employed in the Prandtl wing theory :* 


(;)-% £man(o%) ° 


Substituting in (7) and (2’), we obtain the following expression for the wave drag: 





x pU? N 
R=— I? nb, 9 
-" X (9) 


—again analogous to a well-known expression for the induced drag of a wing.* 
3. Minimum wave drag for given volume and length. The expression for the cross- 
sectional area S corresponding to (8), in the approximation represented by (3), is 


sin (n — 1)0 sin (n + “| 
—1 n+1 f 


It is clear that for closed projectiles, pointed front and rear, b; must vanish.* Also, 
the total volume occupied by the projectile is 


l 1\3 
Vol. -f Sdx = (5) (b; — $62) (11) 
> 


or, for closed pointed bodies, 
l 3 
Vol. = — ins) bs. (12) 


N 


I? 
a {te — 0+ }sin 20), - > | 


2 





(10) 


Hence, for given length and volume, the minimum wave resistance is obtained 
when only J» is different from zero. The geometry of this body is given by 














mlb, L ml*bo . 
S=- oo (sin @ — 3 sin 36) = — sin’ 6 (13) 
and its wave resistance is 
3 U? 
ee et 
| 2 2 
9 U? dmax \* 
otha 
8 2 l 
8 pU?/1\?T Vol. 7? 
2 TaaT : 
e 2\32 (1/2)8 


This is eight times the wave drag of von K4rm4n’s ogive* of equal length and vol- 
ume, or about 11.1 times that of von K4rm4n’s ogive of equal length and caliber. (It 


3 Th. von K4rm4n and J. M. Burgers, Aerodynamic Theory, edited by W.F. Durand, vol. 2, J. Springer, 


Berlin, 1934, pp. 172-175. 

* If b, ~0 while bp =b;= - - - =O, the ogive considered by von K4rmé4n! is obtained. Its maximum 
cross-sectional area is 2*/b,/4 and occurs at its stern, x=/. According to (9), its wave drag is 
(x®/4)(pU2/2)2Pb} or (pU?/2)Smax (dmax/l)*, where Smax is its maximum cross-sectional area and dmax is its 


maximum diameter, or caliber. 
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should be mentioned that this comparison may be misleading in view of the fact that 
von Kaérm4n’s ogive has a blunt stern, so that its wave drag certainly does not repre- 
sent its entire resistance, even in the absence of skin friction. Nevertheless, the wave 
resistance of that ogive may be taken as a convenient reference.) 

The shape of the forward half of the symmetrical projectile represented in (13) 
is drawn in Fig. 1, for the case / =4dmax. For comparison, there is also shown the shape 
of von Karm4n’s ogive having the same caliber and one-half the length. 








Fic. 1. Profiles of various projectiles of minimum wave drag: (a) volume and length given, (b) caliber 
and length given, (c) von Ka4rm4n’s ogive of equal caliber and one-half the length. (Projectiles (a) and 
(b) are symmetrical fore-and-aft.) 


4. Minimum wave drag for given caliber and length. To attack the problem of the 
body shape for minimum wave drag, caliber and length specified, we return to the 
expression for the wave drag given in (2’) and (4) and employ the methods of varia- 
tion calculus. By virtue of the symmetry with respect to x and &, the variation of the 
resistance with varying body form assumes a simple form; viz., 


l l 
6bR = — ro} f ap'(a) f f'(&) log | x — &| dtdx 
0 0 
l l 
+f rf a7@ 10g | x - ¢| dean} 
0 0 
l 
= — 2p f 5f'(x)F(x)dx. (15) 


In this section we shall provide for the possibility [excluded in obtaining (2’’) ] 
that dS/dx, and therefore f(x), is discontinuous at the station where the maximum 
diameter occurs, x =m. Hence, integrating by parts in (15), and again assuming sharp 
points at bow and stern, we write 
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5R = — 2xp f{romal f(x) |m — f af(2)P"(a)ae 
U l 

= — 2xp {Fom6[/(2) In 4 = f i5(a)P"(a)dah, (16) 


where [f(x)]m denotes the value of the discontinuity in f(x) at x=m, and the area 
function S(x) has been assumed to be continuous. 

In the form (16) it is clear that the shape of the part of the body forward of the 
maximum section at x=m can be held fixed while S(x) is varied over the rear part 
to achieve a minimum of R; then the rear shape can be fixed in this minimum-drag 
configuration while S(x) is varied in front to minimize R; the result will be the mini- 
mum-drag shape for given maximum cross section at x =m. We shall also assume that 
the discontinuity of slope represented by [f(x)]m is not varied in the process; it will 
appear later that this is valid. The minimum-drag condition 5R=0 is then obtained 


when 
F’(x) = 0 | 
os s 5 wm, 
F(x) = ox + Ce 
(17) 
F’(x) = 0 
mn S23 il 
F(x) = c3x + 


The analogy with the induced drag of a wing is again useful. The analogous prob- 
lem is the following: to determine the spanwise circulation distribution f(x) so as to 
obtain minimum induced drag, it being required that the total lift be zero, but that 
the lift carried on one side of a station x =m have a given value, equal and opposite to 
that carried on the other side of that station. The result obtained in (17) states that 
the condition of minimum drag results when the downwash F’(x) is constant in each 
of the two parts of the wing. 

Fortunately, investigations have been made* of the behavior of the circulation 
distribution near a point on a lifting line where the downwash in discontinuous. It 
is found that the circulation function is continuous but has a vertical tangent and dis- 
continuous curvature at such a point. Applying this result to our projectile problem, 
we can conclude that f(x) will exhibit a singularity of this type at x =m. Moreover, 
since F(x) can be interpreted as the downwash corresponding to the circulation dis- 
tribution S(x), we conclude that F(x) cannot be discontinuous at x =m if we exclude 
singularities of this type from the shape function S(x). Accordingly, we write 


(c1 ~ C3)m = 4 — Co. (18) 


The resistance of the minimum-wave-drag body is easily calculated from (2’); it is 





‘ A. Betz and E. Petersohn, Zur Theorie der Querruder, Z. angew. Math. Mech. 8, 253-257 (1928); 
also Nat. Advis. Com. for Aeron. Tech. Memo. No. 542 (1929). 

5H. Multhopp, Die Berechnung der Auftriebsverteilung von Tragfligeln, Luftfahrtforschung 15, 153- 
169 (1938). 
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a 
Il 


U 
— xp - (cs — ¢1)Smax + (cxm + 2) [ne] 


T 


r 


U 
p 2 (¢1 ead CS mam (19) 


ll 


where Smax is the cross-sectional area at x =m. 
The form of f(x) corresponding to (17), and subsequently the shape of the body, 
can be determined by inverting (7) by means of the so-called “Reciprocity Theo- 


rem” 


9 


1 2 6-8 
f(x) = —f F(é&) cot = dd. (20) 
47° 6 


The quadratures involved are rather tedious, but can be carried out. The result is 


1 — cos (6 + pz) 
1 — cos (6 — yz) 





f(x) = a5 {( — ¢1)(x — m) log + [(¢s — ¢1)u + re, ]1 sin oh , (21) 
= 

where m = (//2)(1+ cos yw). It can quickly be verified that this function has the type of 

singularity at x =m that was predicted by the wing analogy. 

This expression can be integrated again to evaluate the constants c, and cs; and 
then to determine the function S(x). By integrating from x =0 to x =/ it is determined 
that (¢3—¢:)(u—} sin 24)+7c,=0. Finally, by carrying out the lengthy quadratures 
necessary to apply the condition (U/27)Smax =Jo'f(x)dx, it is found that 


w— $sin 2p 








ty = 4USmax -—— (22) 
I? sin4 u 
The wave resistance (19) then assumes the form 
pU* ~y rr : i 
R= P— Suu( =) ‘ a 
2 l sin4 pu 
or (23) 


pU? Gmax \* m m \?) -? 
ro Palle) bE (BY 
2 l 1/2 1/2 


Thus the wave drag varies symmetrically about m=//2 or n=7/2, and is least if the 
maximum cross section is located at mid-length—i.e., for a symmetrical projectile. 
The wave drag of this projectile is r* times as great as that of von Ka4rm4n’s ogive of 
equal length and caliber. Its shape is indicated in Figure 1. 

5. Concluding remarks. A somewhat similar analysis of projectile shapes for mini- 
mum wave resistance has been made by Haack,’ who considered only symmetrical 
projectiles. The results obtained here are in agreement with Haack’s for such projec- 
tiles, except for the value of the drag of the minimum-wave-drag body for given length 
and volume, which seems to have been tabulated erroneously in the earlier paper. 





6 R. Courant and D. Hilbert, Methoden der mathematischen Phystk, vol. 1, J. Springer, Berlin, 1931, 


p. 83. 
7 W. Haack, Geschossformen kleinsten Wellenwiderstandes, Bericht 139 der Lilienthal-Gesellschaft fiir 


Luftfahrt. 








THE BOUNDARY LAYER IN A CORNER* 


BY 


G. F. CARRIER 
Brown University 


1. Introduction. The laminar flow of a relatively non-viscous fluid through a chan- 
nel is characterized by the presence of a thin boundary layer along the walls. In 
straight channels, such boundary layers are usually assumed to have the velocity 
distribution determined by Blasius [1] for the flow past a flat plate, and the flow 
pattern in the neighborhood of any corner is not mentioned. It seems of interest to 
develop here the change in the Blasius flow implied by such a corner. 

2. The boundary layer problem. We shall consider the laminar flow of an incom- 
pressible fluid which impinges with the uniform velocity V on the edges x =0 of the 
half planes y=0, 2=0. 

The Navier-Stokes equations and the continuity condition which govern such 
flows are 

(v-grad) v + p~! grad p = vAv, (1) 
div v = 0. (2) 


Here v is the velocity with components u, v, w; p is the pressure, v the kinematic 
viscosity, and p the density. 

As v. Karman has pointed out [2], the essence of the treatment of such equations 
in a boundary layer problem is to eliminate higher order terms (by a perturbation 
scheme or otherwise) in such a manner that the order of the equations is not decreased. 
In this way no boundary conditions need be relaxed. We may accomplish this by using 
what is essentially Prandtl’s coordinate transformation [1], namely 


n= y/(vx/V)2, ¢ = &/(vx/V)*!2. (3) 


We also define the parameter £ = (v/ Vx)"/?. 

Since the flow both within and outside the boundary layer may be expected to 
be essentially in the x direction and slowly varying in x, we may attempt to find a solu- 
tion in the form 


u = V[uo(n, £) + Eua(n, $) + Fue +--- | (4) 
v = V(iv, + fe + ---?) (5) 
w= V(éwi + fwe+---) (6) 
b= pVUpot+tpit:--). (7) 


We commence the series for v and w with a term of order £, because we wish a solution 

for which v/ V, w/ V, are small. Furthermore, if we included terms 2, wo, the following 

set of equations would contain terms of order & with no contribution from the 

viscous terms of Eqs. (1) and (2). Thus the solutions wherein vo, wo were not identically 

zero would not provide results corresponding to the phenomenon under investigation. f 
* Received Aug. 30, 1946. 


t Actually, the fact that our results constitute a solution which obeys the differential equation and 
boundary conditions is sufficient justification for taking 19 =w»e=0. 
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The substitution of Eqs. (4) to (7) into Eqs. (1) and (2) leads to the system 





te) te) 
-= (nduo/On + £0u0/OC) + 010U0/dn + w1Ou0/dB — = = - 2 = 
0? 0? 
— (Sot sa)m tHe 20 (8) 
= AW fe e(S| ss ~"|+ tn )+ =0 (9) 
an an ct im. & On? 7 
Opo Opi 
x + a £( + (10) 
Ou ¢ Ou dv Ow 
jp PE a a a En Pr 


The solution of this system of equations requires that the coefficient of each power 
of £ in each equation vanish. The first order approximation to the result is defined by 
the vanishing of the coefficients of £°. The result can be expected to be valid only when 
the remaining terms of the series are negligible, that is when & is small. Thus the 
solution, like that for the flat plate, is valid only at sufficiently large distances from 
the leading edges of the planes. 


We now note that the £° terms of Eqs. (9) and (10) vanish only if »=const; the 
£° term of Eq. (11) vanishes if we write 
uo = Sar(n, $), = 1 = F(nBax — Br), =—§- Wi = «HSB — 8 a)- 
Thus it remains to find g(n, ¢) such that, 
g(0, ¢) = g,(0, ¢) = g(n, 0) = gr(n, 0) = O 


and 
lim Sar(n, o) eas i, 


9 $00 
the implied symmetry condition 
g(a, b) = gb, a), 
and the differential equation implied by Eq. (8) 
Saat + Sartre + Fl erSonr + Sagar} = 0. (12) 


We may expect that far from the corner the solution will be essentially that for 
the flat plate. Hence, we write 


g(a, £) = fon) fo(S) + h(n, $) (13) 


where fp is that solution of 


2f"" + ff’ = 0 


such that f(0) =f’(0) =0; fd...(a@) =1. This function is tabulated in [1 ]. 
Equations (12) and (13) lead to the equation 
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honnt + Angee + 20(n, $) Agee + 20(8, 2) Aare 
+ b(n, §) h,+ b(f, n)hy + $(hghyne . hyhger) = asA(n, $), (14) 


where 
a(n, £) =dfoln) fol), b(n, $) = fo’ (nfo (S) 
A(n, £) = Ef folm) fe’ (mf (O[1 — 6] + fol fe’ Ose (n) [1 — fe (a) ]}. 


This equation may be integrated once each over 7 and ¢ taking account of the 
boundary conditions to yield (when g¢= —25h,;) 


1 t 
Ag + 2a(n, £)dg/dn + 2a(f, n)de/dF + b(n, of edn + b(f, vf edf 
0 0 


1 . t 
tole gdn + ef oat | = gs A(n, £). (15) 


The boundary conditions are 


(0, g) ad e(n, 0) = lim ¢(n, §) = 0. 
n $70 
This last form of the equation seems best suited for numerical evaluation. The relaxa- 
tion method [3] appears to be the most appropriate for the determination of ¢ so 
we form the difference equation derived from Eq. (15) by taking points spaced unity 
apart in 9 and ¢. The subscripts m and n are used to index these point positions. The 
difference equation is 


Pm+1,n + Pm—1,n + Pm,n—1 + Om,n+1 — 4Qmn + Sudl@iust.o — Om—1,n) 
n m 


+ Canl On wed — a + ban f gdn + Pie gdt 


0 0 


+ .01 | (omc * emt) f edt + (Gm+i,n = emi) f edn | + Awa = 0. (16) 
0 0 

In this equation the integrals may be evaluated by the simple trapezoidal rule since 

the function ¢ is very “smooth” although if more accuracy is desired a simple graphi- 

cal method is conveniently employed. 








TABLE I 
¢ | $(n, £) 

rn | 0 1 2 3 4 5 6 7 8 fo (n) 
0 | oO 0 0 0 0 0 0 0 0 0 

1 | 0 58 1.00 1.00 64 .25 08 .02 .00 .330 
>i «s 1.00 1.60 1.46 .86 .28 .08 .02 .00 .630 
3 | 0 1.00 1.46 1.23 61 .16 .04 .01 .00 846 
4 0 .64 .86 61 .24 .03 .01 .00 .955 
5 0 .25 .29 .16 .03 01 .00 .992 
6 0 .08 .08 04 01 .00 00 999 
7 0 .02 

~ 0 .00 
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The numerical procedure is this: guess values for ¢ at all points m, n=8. Replace 
the zero on the right side of Eq. (12) by Qn, and compute each Q,,, (the residuals). 
Then revise the guesses for the gm» in such a way as to decrease the Qn», disregarding 
the changes in the values of the terms containing integrals. When considerable im- 
provement has been made, recompute the Q,,, using the complete equation (12) and 
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Fic. 1. Contours of constant U in corner boundary layer. 


repeat the foregoing procedure. It is not necessary to get extremely accurate values 
of g (especially since a, b, A are not known too finely) because the velocity uo =f¢ (7) 
fo ()+es¢(n, £) will be accurate to three places when ¢ is known to the one hun- 
dredths digit. The functions ff and ¢ are tabulated in Table I and contours of con- 
stant uo are shown in Fig. 1. 
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THE TREATMENT OF SINGULARITIES OF PARTIAL 
DIFFERENTIAL EQUATIONS BY RELAXATION 
METHODS* 


BY 


H. MOTZ 
University of Sheffield, England 


1. Introduction. In the course of a study of boundary value problems arising in 
radiation theory and electrostatics, the treatment of singularities demanded special 
attention. In most problems of practical importance boundaries with sharp corners 
occur. Such sharp corners give rise to singularities of various types. When the com- 
puted function is bounded, but has a branch point at the corner, the difficulty is 
not serious. The use of a graded net with a finer mesh size near the corner is possible. 
Conformal transformation which automatically provides a finer net near corners is 
also successful. The mesh size near the corner should be of the order of magnitude of 
the radius of curvature of the corner, and when this is small a mathematical idealiza- 
tion involving infinitely sharp corners is preferable. The special treatment outlined 
in this note makes use of such an idealization and shortens the labour considerably. 
Special treatment is essential when the function approaches infinite values near the 
corner. 

2. Plane harmonic functions. Solutions of ¥26 =0 are bounded when the boundary 
condition prescribes constant values near the corners. It can be shown that they are 
also bounded when 0¢/0»v is constant, where v is the direction normal to the boundary. 
This type of boundary condition occurs e.g. when two plane harmonic functions @ 
and W are computed inside a boundary B for the purpose of a conformal transforma- 
tion 

x+iy= $(x, y) + ip(x, y) 
and yY=const. is specified at the boundary forming the corner. When expressed in 
polar coordinates r, 3 centered at P (Fig. 1), the equation 





Vv’ = 0 (1) 
becomes 
0? 0 0? 
in bir henet (2) 
or? Or =a? 
With (r, 3) =R(r)- O(8), the following equations for R and @ are obtained 
ao 
+ nO = 0, (3) 
dv? 
1?R dR 
Po gip lane ye ow & (4) 
dr? dr 


In these equations m? stands for a positive constant, and 


dO 
—=0 when 8=0, J=a. 
dd 


* Received Jan. 2, 1946. 
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Hence . 
¢= ‘» Axr" cos nd (5) 


kem—co 


where n=mk/a (R=0, +1, +2, +3,---). 
In order to investigate the terms with negative exponent in this series, we exclude 
the corner by a small circle of radius p. On this circle (0¢/07r),.0=0. It is found that 


A..= p44, (s=1,2,3,°-:). 
When p-— 0, the circle contracts towards the point P and the terms with negative 


exponents vanish. Thus ¢ will be represented by the series 


T 2x 3nr ps 
@ = Ao + A;r™!* cos —B + Aogr?*/* cos —d + Azr**/* cos—B+:-- (6) 
Qa 


a a 


in the neighbourhood of P. 





Fic. 1 


3. Method of special treatment. The method of treatment will now be explained 
with reference to the example of a corner with a=27, r/a=}. In the treatment of 
two-dimensional problems by relaxation methods,’:? the function ¢ is computed at 
points of a net with small but finite mesh size. Let us denote by ¢o the value of @ 
at such a point, by gu, ¢2, $3, ds the values of ¢ at the nearest neighbouring points. 
The mesh length is a. At points where the function is regular, double Taylor expansion 


shows that 


4 1 a 1 ae 
a’V*d — z dm + 460 = — D Ao (x) - DR Ao(y)---, (7) 


m=1 


where Ai¥(x), Ai’(y), are the fourth central differences in the x and y direction, re- 
spectively, at the point where ¢ =¢o. This expansion can only be used when the right 
hand side converges. At a singularity and its nearest neighbouring points this expan- 
sion is not valid. Figure 2 shows an example of a boundary where ¢=0 on AE, 
¢ = 1000 on EB, and 0¢/dv =0 on all other boundaries. The Taylor expansion fails at 





1H. W. Emmons, Numerical solution of partial differential equations, Quarterly Appl. Math., 2, 


173-195 (1944). 
2 D. N. de G. Allen, D. G. Christopherson, L. Fox, J. R. Green, H. Motz, F. S. Shaw and R. V. South- 
well, Relaxation methods applied to engineering problems, Phil. Trans. Royal Soc. London (A), 239, 367- 


386, 419-537, 539-578 (1945). 
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P’, QO’, R’, and S’. In order to obtain valid equations at these points, we consider series 
of the type (6) at the pivotal points P, Q, R, S 


v 
& = Ao + Ajr'!/? cos - + Aor cos 3 + Azr?/? cos 3/28 +-°:-. 


(8) 


Only the first four terms are retained. The units of r can be so chosen that r=1 at 


P, O, R, S. In terms of $p, de, dr, Os one obtains 





















































Cc 
589 598 625 671 737 818 907 iddo 
581 58 615 661 730 814 905 1004 
P & 
557 563 585 6 708 803 901 1004 
+10 +13 
Pp! Q' 
520 523 533) 561 669 789 89% 1004 
+14 
F 
474 475 464 436 328 208 10 0 
+1 Rt s! -14 
44d 434 412 368 £89 194 97 ) 
R -10 -13 3 
41 40 38 336 267 183 93 0 
40 39 37 326 260 17 9 0 
D A 
Fic. 2 


0.25 (6p + de + or + $s), 
0.191(¢p — dr) — 0.462(ds — $¢@), 
0.354(¢g — op — or + $s), 

— 0.191(¢e@ — os) + 0.462(¢r — ¢$p). 


At the special points P’, Q’, R’, and S’, we find from (8) 
0.457p + 0.23549 + 0.2096r + 0.09995, 
0.235¢p + 0.593¢eq + 0.099¢r + 0.073¢s, 
0.2094p + 0.099¢9 + 0.4576 + 0.23545, 





op: 
oq’ 
or’ 
os’ = 


Il 


ll 





0.0994p + 0.073dq + 0.235¢r + 0.593¢s. 


(9) 


(10) 













374 H. MOTZ [Vol. IV, No. 4 


These are the equations used at special points. The relaxation procedure is carried out 
normally everywhere, observing that equations (10) hold at special points. The resid- 
ual at a special point due to an increment at a pivotal point P is therefore the prod- 
uct of this increment with the coefficient of @p in the equation which holds at the 
special point. This is in accordance with the usual relaxation procedure. The removal 
of a residual at a special point is particularly easy. It is simply subtracted from the 
value of @ at the special point in question. Due to this removal the usual residuals 
accrue at ordinary neighbouring points, but of course, no residuals are passed on to 
special points. 






































YG 
Li 


Figure 3 refers to an example with a=3-2/2. Here we retain three terms only 
Special points are P’, Q’, R’, pivotal points P, Q, R. The equations at special points are 


op: = 0.486¢p + 0.257¢9 + 0.2575, 
oq = 0.2576p + 0.61269 + 0.131¢s, (11) 
dr’ = 0.257op + 0.131¢¢ + 0.612¢s. 


It should be checked whether the three first terms of the series 
@ = Ao + Ayr*!® cos $0 + Ar24!* cos 40 + --: (12) 


which has been used for the derivation of (11) represent the function ¢ adequately. 
This is done by comparing the result of the relaxation computation at points S, 7, 
where ordinary difference equations have been used with the values of ¢ calculated 
by means of the first three terms of (12). 

A similar check was carried out at analogous points in the example of Fig. 2. It 
was found that the agreement was not satisfactory. The errors have been recorded in 
Fig. 2 underneath the respective @ values. In this case it is possible (with the net 
shown in Fig. 4) to retain five terms of the series. Pivotal points are 7, U, V, W, X; 
special points 7’, U’, V’, W’, X’, and the equations for ¢ at special points are 


gr’ = 0.546 dr + 0.313 du + 0.062s¢y + 0.0625¢w + 0.016¢x, 
gu: = 0.156 or + 0.578 du + 0.188 dy + 0.047 dw + 0.031 dx, 
gv’ = 0.031 dr + 0.188 du + 0.562 dy + 0.188 dw + 0.031 dx, (13) 
dw: = 0.031 dr + 0.047 du + 0.188 dy + 0.578 ow + 0.156 dx, 


0.016 dr + 0.062s6u + 0.062 dy + 0.313 dw + 0.546 x. 


ox’ 
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The result is again checked by comparing the result of the relaxation procedure with 
the @ values near the corner calculated from 


= Ao + Ayr’? cos 38 + Aor cos & + Azr*/? cos 48 + Agr? cos 28, 

































































where the A’s are given by 


A, 


0.250(¢u + dv + ow) + 0.125(¢r + $x), 
0.354(¢uv — ow) + 0.250(¢r — x), 

— 0.500¢y + 0.250(¢r + ¢x), 
0.354(¢w — ou) + 0.250(¢r — ox), 
0.250(¢v — du — dw) + 0.125(¢r + ox). 


The agreement is now much better. In Fig. 4 the errors have been recorded. It will 
be noticed that the mesh points of Fig. 2 lie between those of Fig. 4. By interpolating 
values at midpoints of the meshes of Fig. 4, we find that the solutions given in the 


two figures are in fair agreement. 


When the above test fails, a finer net should be used as a rule, because the calcula- 
tion becomes rather cumbersome when more than five terms of the series are retained. 
To obtain, without the special treatment, a result which differs from the one of 
Fig. 2 by less than 1% at any point of the net, the net near the corner would have to 


be 7 times as fine. 


590 607 643 700 776 863 954] —flooo 
U 
573 589 623 682 765 857 952] 000 
d : 
641 55 578 639] 745 849 950| L000 
atl al 
v vt tT T 
9 49 499 499] 84 949] _ ooo 
271 156 51] |o 
x! x 
wi. 
457 446, “ 359 25 150) 50) =p 
+ + 
W 
45 409 375 316 234 1 aa so 
396 35 137 0 
Fic. 4 
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(14) 
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4. Other examples. As an example of a corner where the value of the function is 
specified, let us consider an electrostatic potential y. In this case the series is 


YW = Ao + A;r/? sin 48 + Ag sin 3 + Ar®/? sin FV + --- | (15) 


The components of the electric field in Cartesian cordinates, E,, Ey, are not bounded 
on a sharp corner when a>z. Let us consider the term r* cos nd of the series (5). 
E,and E, will contain terms r*~ sin (n — 1), r"—! cos (n—1)8, respectively, and nega- 
tive exponents of r will therefore occur when a>. The method outlined above can 
still be used to compute a function with such a singularity. In the case a=27 the 
negative exponent —} occurs. Terms with exponents —} and +} depend on @ in the 
same manner. It is therefore necessary to choose pivotal points which have not all 
the same distance from the corner. 
The method is equally applicable to solutions of the wave equation 
V*o + k*o = 0. (16) 

In Cartesian coordinates and with § =x/a, 7 =y/a (where a is the mesh length of the 
net), Eq. (16) becomes 

a*p da 

— + — + k*a% = 0. (17) 

ag an? 


Referring again to the case a =27, 0¢/dv =0 and retaining five terms, we see that the 
expression (5) holds at pivotal points. At the special points we have 





Jio( Ra J (ka J 3/2( ka) 
Pee ee ee ee Pena dail tdi tea te 
J 12(2 ka) J ,(2ka) J 3;2(2 ka) 
Jo(k 
4 Aan cos 28, 
J2(2ka) 


where J, are the Bessel functions of order n. When ka <0-1, the ratios J,(ka)/J,(2ka) 
differ from (})” only in the third decimal. When the mesh length a is small compared 
with the wave length /=27/k, the special equations are therefore the same for solu- 
tions of the wave equation and those of Laplace’s equation. 

5. Conformal transformation. When a solution of more complicated differential 
equations, e.g. the equations of viscous flow, or V‘F =0, is computed it is often an ad- 
vantage to remove singularities at the boundary by a conformal transformation 
¢=(x, vy), Y=W(x, y). Let us suppose that it is desired to transform the interior of 
the boundary shtwn in Fig. 2 into the interior of a rectangle in the ¢, y plane. 

The lines @ =const. at suitable intervals can be found from the ¢-values recorded 
in Fig. 2. The lines y’ =const. are orthogonal to the lines ¢ =const. and are best com- 
puted separately. 

The condition for orthogonality 


Op Oy’ dd dy’ 
Ox Ox dy dy 


is satished when 


ad 
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0g oy’ Op oy’ 
—=)—, — = -\— (18) 
Ox oy Oy Ox 


where A is a constant. From these equations it follows that 


2 2 2p’ 4 
-- “~~ = 0, 7. + 7 = (). (19) 
a éy® Ox? dy? 
The boundary conditions for ¥’ are dy’/dv=0 on EA, EB, py’ =0 on EF, yp =const. 
on AD, DC, and CB. The last constant is arbitrary and may be given a convenient 
value, e.g. 1000 for three figure accuracy. 
It is easily seen that it is necessary to determine the constant \ in (18), in order 
to carry on with the computation of the original equation (e.g. V‘F =0). In the co- 
ordinates @ and Y=hy’, this equation becomes 


3 3 o? 
(= + -,) (= + =F = 0. 
dg? ay?] \ag? oy? 


The constant \ is determined by (18). These equations can be regarded as an esti- 
mate of \ at every point. Denoting the finite differences in the x and y directions at 
the mesh point 7 by 


Di), Dwi), Dwi), Dw'(i), 
the quantities 6,(7) and 62(z) defined by 
b(t) = D.p(t) — DD W'(t), a(t) = Dyo(i) + ADW'(2) (20) 
are not all zero, but constitute a measure of the computational error. It is desired to 


find a mean value for \ for which the variance of the computational error is a mini- 
mum. The sum 


2.. 2, 
DX [61(é) + 62(4)] 
is thus minimized with respect to \ and the following expression for \ is obtained: 


> {D.4()DW(i) — Dw()Dw'(i)} 


DX ({DwH (i }* + [Dwi }*) 


It has been found that, with the help of this technique of separate computation of 
the two transformation functions and using the special treatment of corners at the 
boundary conformal transformations can be computed with great accuracy. 

6. Acknowledgment. The numerical work for this paper was done by Miss 
L. Klanfer whose services were put at the author’s disposal by the Directorate of 
Scientific Research of the British Admiralty. 





(21) 











378 


THE GENERAL VARIATIONAL PRINCIPLE OF THE THEORY 
OF STRUCTURAL STABILITY* 


BY 


WILLIAM PRAGER 
Brown University 


1. Introduction. This paper is concerned with the general problem of structural 
stability in the elastic or plastic range. Two slightly different formulations of this 
problem are found in the literature. According to the first, one considers a deformable 
body which, initially, is free from stresses, and which is then subjected to a system 
of loads of gradually increasing intensity. As long as these loads are sufficiently small 
the equilibrium configuration which the body assumes under their influence will be 
stable; one asks for that intensity of the loads for which this equilibrium configuration 
first becomes unstable. According to the second formulation of the problem of struc- 
tural stability, one considers a given configuration of a deformable body and an equi- 
librium system of body and surface stresses and asks whether, in the presence of these 
initial stresses, the given configuration is stable or not. This second point of view is 
adopted in this paper because: 

(1) it clearly separates the stability problem from the problem of finding the 
stresses produced by the given loads, and 

(2) the manner in which the initial stresses are produced is irrelevant for the solu- 
tion of the stability problem. In particular, it is by no means necessary that the initial 
stresses are produced by loads which are applied to an otherwise stressfree body; 
they may be produced by temperature changes or may partly be due to previous over- 
straining of the body. 

-Once this second point of view is adopted, stress-strain relations enter into the dis- 
cussion at one point only: we must be able to predict the infinitesimal changes in 
stress which correspond to the infinitesimal strains associated with a system of in- 
finitesimal displacements from the considered equilibrium configuration. As the rela- 
tions between these infinitesimal changes in the stresses and strains are essentially 
linear, the only difference between the elastic and plastic ranges consists in the fact 
that in the plastic range a different set of coefficients must be used in these linear 
relations according to whether the change of stress constitutes “loading” or “unload- 
ing,” while no such distinction need be made in the elastic range. 

In Section 2, the general problem of structural stability is reduced to an eigen- 
value problem for the displacements from a configuration of indifferent equilibrium 
to a neighbouring configuration of this type. Except for the consideration of plastic 
deformations, we follow Biezeno and Hencky’ in this derivation, but simplify the dis- 
cussion by the systematic use of tensors. In Section 3, a variational principle is derived 
which is equivalent to the eigen-value problem formulated in Section 2. As an ex- 
ample for the application of this principle, the lateral buckling of an unevenly heated 
lamina is treated in Section 4. 


* Received June 5, 1946. 
1C, B. Biezeno and H. Hencky, Proc. Roy. Acad. Amsterdam, 31, 569-592 (1928). 
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2. The eigen-value problem associated with the general problem of structural sta- 
bility. We consider a given configuration of a deformable body and an equilibrium system 
of body and surface stresses which is given to within an arbitrary factor \. Jf Xd is 
sufficiently small, this equilibrium configuration will be stable; we ask for that value of 
for which it becomes indifferent, assuming that the additional stresses which are produced 
by infinitesimal displacements from the given equilibrium configuration are linearly re- 
lated to the corresponding infinitesimal strains. This critical value of X will be called 
the safety factor of the considered equilibrium configuration. With respect to a system 
of rectangular Cartesian coordinates x;, let us denote the components of the given 
stresses by \o;; and the components of an infinitesimal displacement from the given 
equilibrium configuration by u,. If the unit vector along the outward normal to the 
surface is denoted by nj, the surface stresses are 


AT; = doi sni. (1) 

The quantities ¢;; must satisfy the equilibrium conditions 
oi3, = 0, (2) 
where the subscript i after the comma denotes differentiation with respect to x;, and 


the usual summation convention regarding repeated subscripts is adopted. 
The infinitesimal strain associated with the displacements u; is given by 


€¢j = 3(ui,; + U;,). (3) 


Since the relation between this strain and the corresponding additional stress 7;; is 
assumed to be linear, we have 

735 = Cisjnr€nt, (4) 
where C;;x; is a fourth order tensor which is symmetric with respect to 7 and j and 
with respect to k and /. If, in particular, 7;; and €;; are assumed to be related to each 
other by the generalized law of Hooke, we have 





Cijxi = 26¢{ udu _ isbn), (Sa) 


— a 
where Go denotes the modulus of rigidity, v Poisson's ratio, and 6;; is the Kronecker 
delta. If the body under consideration can be expected to behave like an isotropic 
elastic solid for an infinitesimal displacement from the given equilibrium configura- 
tion,” i.e. if the stresses \o;; do nowhere exceed the elastic limit of the material, the 
expression (Sa) may be used in connection with the stress-strain relation (4). On the 
other hand, where the stresses \o;; exceed the elastic limit, different expressions must 
be used for C;;,: according to whether the stresses 7;; associated with the strains €;; 
constitute “loading” or “unloading.” We reserve the complete discussion of suitable 
stress-strain relations beyond the elastic limit for another paper and give but one ex- 
ample here. Defining the stress deviation as sij=0i;—40%%5;; and its intensity as 
S = 4545553, we set 


2M. A. Biot [J. Appl. Phys., 10, 860-864 (1939) ] and, more recently, F. D. Murnaghan [Proc. Nat. 
Acad. Sci., 30, 244-247 (1944) ] have pointed out that an elastic solid under initial stress can be strictly 
isotropic only if the initial stress is of the nature of a hydrostatic pressure. For the conventional structural 
materials, however, this small anisotropy caused by the initial stress can be disregarded as long as the 
initial stress does not exceed the elastic limit. 
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Vv Go = G 
Cijet = 2Gof 5545 72, — ——— 56k} — ——— 5:,Se2 for 5,3¢;; > 0 (Sb) 
1 — 2p S 
and 
Ci jet = 260 548; = iu) for Si j€ij < 0. (5c) 
— 2 


Here Go denotes the value which the modulus of rigidity assumes in the elastic range, 
while G=G(S) is the so-called tangent modulus of rigidity. In the elastic range G=Go, 
and (5b) as well as (5c) reduce to (5a). The stress-strain relations which are obtained 
by substituting (5b) and (Sc) into (4) were suggested by J. H. Laning in an unpub- 
lished paper (1942); they constitute a generalization of stress-strain relations which 
the present author had used in earlier papers.* We note that Cijx1=Cxii;, according 
to (5a), (5b), and (5c). 

A generic particle with the coordinates x, in the initial state has the coordinates 
&;=x,+4u; in the considered neighbouring state, and 


dk; = (655 + uy, s)dx; = (655 + 65 + w;j)dx;, (6) 
where the deformation €;; is defined by (3) and 
Wis = 3(Mi,5 — U;,s) (7) 
is the rotation associated with the displacement 4. 
The infinitesimal force Adf; which is transmitted across the surface element dS in 
the initial state equals 
Adf; = AT;dS = ho; nd. (8) 
The force which is transmitted to the corresponding material element in the neigh- 
bouring state will be written in the form 


df; = Ao; jnidS. (9) 


Note that the normal vector m; and the area dS in the initial state are used in (9). 
This means that the stress tensor \@;; is defined in the Lagrangian manner‘ with the 
initial state as the state of reference. Consequently, A¢;; is not a symmetric tensor; 
it will be written in the form 


- , 4? 
AGizg = AOEs H Tig + Tig H ij (10) 
where the terms 7;;, tj, and 7 are infinitesimal changes of stress defined in the 


following manner: 

(1) the tensor 7; is symmetric; it represents the change of stress associated with 
the infinitesimal strain €;; and is given by Eq. (4); 

(2) the tensor 7, too, depends on the strain €;;; it is antisymmetric and repre- 
sents the change of stress necessary to restore the moment equilibrium which is 
expressed by the symmetry of o;; in the initial state and which is disturbed by the 


deformation; 





3 W. Prager, Proc. 5th Internat. Congr. Appl. Mech. Cambridge, Mass., 1938, pp. 234-237; Prik- 
ladnaia Matematika i Mekhanika 5, 419-430 (1941); Duke Math. J. 9, 228-233 (1942). 

‘H. Jeffreys has recently given a similar analysis using the Eulerian approach [Proc. Cambridge 
Phil. Soc. 38, 125-128 (1942) ]. The Lagrangian approach seems more suitable, however, for the problem 
under consideration. 
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(3) the term 7;}, finally, depends on the rotation w,;; it represents the change of 
stress, with respect to the fixed coordinate axes, which is produced by this rotation. 

Since only first order terms in €;; and w,; need be considered in the following analy- 
sis, the order in which the deformation €;; and the rotation w;; are applied is imma- 
terial. 

The antisymmetric tensor 7j, depends only on ¢€;;. To find its mathematical ex- 
pression, it is therefore sufficient to consider a pure homogeneous deformation, i.e., a 
deformation for which u;,; is independent of the coordinates and u;,;=4;,;= €4;. On 
account of (9); the equations of equilibrium for the deformed body are 


J aanas = 0, feast - Cink; )ndS = 0 


f testo = 0, ff Gsm _ Gi-% 5) dv = (). 


Since these equations must hold not only for the entire body, but also for an arbitrary 
portion of it, we must have 


or 


O85. = 0, (11) (65 ;%% Zs OsrXj),s = 0. (12) 
For the considered pure deformation, 7j;=0 and 
Hi,5 = O55 + us,5 = O55 + sj. 


Using the symmetry of the tensors o,; and 7;; in addition to the Eqs. (10), (11), (2), 
and neglecting higher order terms, we may therefore write (12) in the form 


Thi 7 Thi = 2ri; = Nosnens = O jk€ki)- (13) 


The tensor 7{; depends only on w,;. To find its mathematical expression, it is suffi- 
cient to consider a rigid body rotation, i.e., a system of displacements u; which de- 
pend linearly on the coordinates x; and satisfy u;,;= —u;,;=@i;. By this rotation the 
components of the infinitesimal force transmitted across a given surface element are 


transformed according to 
df; = (555 + ui,s)df; = (655 + wis dfs = dfs + wisA fi. (14) 


For the considered rigid body rotation 7;;=7;;=0. Using (8), (9), and (10), we 
may therefore write (14) in the form 
vt (15) 


| alae AG KWE j- 


Returning now to the consideration of arbitrary infinitesimal displacements 1, 
we write in accordance with (10), (13), and (15): 


Giz = Aoiz + rig + FACOG — 7 jRERE) — ATERWR; (16) 
On account of (2), the equilibrium condition (11) furnishes therefore 
[rij + §A(osnens — o jens) — Aowes],s = 9, (17) 
and the condition df;=df; furnishes 


[rs5 -f EA(OsKER; —_ O jx€ki) - Aosnwn 5] = 0. (18) 
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Except for our more general definition of the tensor 7;;, Eqs. (17) and (18) agree with 
those derived by Biezeno and Hencky. Biot® obtained the same relations from his 
non-linear theory of elasticity, and Neuber® has recently discussed the formal rela- 
tion of the differential equations (17) to the fundamental equations of elasticity. As 
was already pointed out by Biot, Eqs. (17) differ somewhat from the equations which 
Trefftz’? derived using an unconventional definition of stress. If the given state of 
stress, \o;;, is homogeneous and if the coordinate axes have the directions of the prin- 
cipal axes of this state of stress, Eqs. (17) reduce to the form given by Southwell.’ 

By means of (3), (4), and (7), the quantities €,;, 7;;, and w;; can be expressed in 
terms of the first derivatives of the displacement u;. In this manner an eigen-value 
problem for the displacement u; is obtained. The smallest eigen-value \ is the desired 
safety factor for the given distribution of initial stresses. We refrain from formulating 
this eigen-value problem explicitly, because in all but the most simple cases its exact 
solution would hardly seem possible. 

3. The variational principle associated with the general problem of structural sta- 
bility. The form of Eqs. (17) and (18) suggests the existence of an equivalent varia- 
tional principle from which approximate solutions of stability problems can be ob- 
tained. Indeed, let us establish the Euler equations and natural boundary conditions 


of the variational problem 
sf [Cpers€pgtre + AT pq( Ur, pltr.g — Erp€rg) ]dv = 0, (19) 


where only the displacements u, and hence strains €,, are to be varied, but not the 
stresses ¢,, and the coefficients C,,,, which depend on the stresses. If the integrand 
of the left-hand side of (19) is denoted by F, the Euler equations and natural bound- 


ary conditions are 











a / OF oF 
paste = 0, (20) n;, = 0. (21) 
dx;\0u;,i OU 5,4 
Since 
O€ng 
ee 3 (5 jig + 5p) ja), 
OU j,4 
we have 
OF . 
= 2C 3 je1€Kt + A[2osuu;,1 — O¢h€ jk — o 5jx€it | 
OU ;,: 


- 2r5j + A [oie ; —- €¢ ain 265404 5]. 


Equations (20) and (21) thus are indeed identical with (17) and (18). 

The variational principle (19) can be used in very much the same manner in which 
the principles of minimum potential energy and minimum complementary energy 
are used in elasticity:* by reasonable assumptions concerning the displacements 1; 





5M. A. Biot, Phil. Mag. (7), 27, 468-489 (1939). 

6H. Neuber, Z. angew. Math. Mech. 23, 321-330 (1943). The author is indebted to Professor 
E. Reissner for the reference to this paper. 

7E. Trefftz, Z. angew. Math. Mech. 13, 160-165 (1933). 

* R. V. Southwell, Phil. Trans. Roy. Soc. London (A), 213, 187—244 (1913). 

® See, for instance, E. Volterra, Atti Accad. Lincei, Rend. (6), 20, 424-428, 463-467 (1934); 21, 14-19 


(1935); 23, 329-332 (1936). 
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the class of admitted functions is restricted and the variational problem simplified. 
In using this technique, we must see to it that the restrictions imposed on the displace- 
ments u; do not rule out the possibility of fulfilling the boundary conditions (17). 

4. An example. To illustrate the manner of application of the variational principle 
formulated in Section 3, let us discuss the lateral buckling of an elastic, prismatic 
beam of the length / which is built in at both ends. We assume that the cross section 
of this beam is doubly symmetric. Taking the origin of the coordinates at one end 
of the beam, we let the axis of x; coincide with the axis of the beam and the axes of 
x2 and x3 with the axes of symmetry of the cross section x;=0. To simplify the expres- 
sion (5a) for the coefficients Cij;x1, we shall assume that y=0. This assumption is in 
conformity with the spirit of the engineering theory of the bending of beams; in 
using it we must keep in mind that Young’s modulus Ep» equals twice the modulus of 
rigidity Go if v=0. 

As to the initial state of stress, let us consider the case where 


O11 = CX2, (22) 


while all other components of o;; vanish. The constant c in (22) obviously has the 
dimension of a stress divided by a length. In an originally unstressed beam with 
built-in ends a stress distribution of the type (22) can be produced by changes of 
temperature which vary linearly with x2. If the width of the beam (measured in the 
direction of x3) is small in comparison to its height, (measured in the direction of xs) 
the stresses (22) may produce lateral buckling. The infinitesimal displacements as- 
sociated with this type of instability may be described in the following manner: 
a generic cross section x; of the beam undergoes a translation “(x;) in the direction 
of the x;-axis, a rotation —u’(x,;) about the x-axis which makes the cross section 
remain normal to the bent centerline of the beam, and, simultaneously, a rotation 
—6(x;) about the x-axis; in addition to this rigid body displacement the cross section 
undergoes a warping —w(xe, x3)0’(x1) which is associated with the twist —6’(x:). The 
corresponding displacement components are 


My = — xX3u'(%1) — w(x2, %3)0’(x1), ue = X30(2), us = (x1) — x2(%1). (23) 


Note that on account of the assumption vy=0 the longitudinal extension 0u/0x, is 
not accompanied by any lateral contraction. Particularly simple expressions for uz 
and wu; are thus obtained. The matrices of the derivatives u;,; and of the strains ¢€;; 
therefore are 








’ — x3u’’ — w6”’ — 0 dw/dxs — u' — 0 dw/dx; 7 

“ui; = x30’ 0 —@0 , (24) 
x u! — x26" 0 0 
- — x3!” — wh" 10’(x3 — Ow/ON2) — 36'(x2 + Ow/dx3) 7 

ij; = 40'(x3 — Ow/dxe2) 0 0 ; (25) 
L— 36’(x2 + Ow/dx3) 0 0 od 


Since o;;=0 unless i=j=1, we need only uj —€xn€x for the evaluation of the term 
with the factor \ in (19). Now, for a doubly symmetric cross section the warping func- 
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tion w is odd in x2 as well as in x3. Taking account of this fact, and keeping in mind 
that oy is odd in x2 and even in x3, we find that 


l 
fesse — €41€4;)dv = — 2c f wet side = — dels f u’0'dxy, (26) 
0 


where J; denotes the moment of inertia of the cross section with respect to the x3-axis. 
We now proceed to the evaluation of term Cyors€pg€re in (19). With »=0, Eq. (5a) 

takes the form Cijx2.=2Go06:,6;, and the stress-strain relation (4) reduces to 
Ti 2G oe; j. (27) 


In applying this, we shall replace 2Gp by Eo whenever 1=7. In view of (25), we have 


Cpars€patre = Tpo€pq = Eo(xgu’’ + w0’’)* + 4Go(ej2 + €j3), (28) 


where €,. and €3 depend on the twist 6’ and on the warping w per unit twist in pre- 
cisely the same manner as in the case of pure torsion. In this case, however, the in- 
tegral of 4Go(€.+ 3) over the cross section equals GoC8’2, where GoC denotes the tor- 
sional stiffness of the beam. Adopting the warping w per unit twist found in the case 


of pure torsion, and setting* 
l= f wd A, (29) 


where dA denotes the area element of the cross section, we obtain 


l 


l l 
f Cpore€poérsd? = Eol f ul*dax, + El f 6’dx; + GoC f 6’2dx1, (30) 
0 0 0 


where J, is the moment of inertia of the cross section with respect to the x-axis. 
Substituting the expressions (26) and (30) into (19), we obtain 


Eoloul” + rcI30” = 0, EO — G.CO” + rcIyu” = 0 (31) 


as the Euler equations for our problem, and 
6’ =(Q for 4=0 and 4 =/ (32) 


as the natural boundary conditions. In addition to these natural boundary conditions, 
we have the imposed boundary conditions 


dG=u=n’'=0 at x4,=0 and x, =i. (33) 


The safety factor \ is found as the lowest eigen-value of the problem formulated by 
Eqs. (31), (32) and (33). 


* Note that for the doubly symmetric section considered here the point x;, 0, 0 is the shear center 
of the cross section x;. Since w is odd with respect to x, and x3, we have w=0 at this point. These remarks 
identify the definition (29) with that given by J. N. Goodier, Eng. Exp. Station, Cornell University, Bulle- 
tin No. 27 (1941), p. 9. 
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UNSTABLE SOLUTIONS OF A CLASS OF HILL 
DIFFERENTIAL EQUATIONS* 


BY 


GABRIEL HORVAY 
McDonnell Aircraft Corporation 


1. Introduction. Linear differential equations with periodic coefficients play an im- 
portant role in problems of engineering and physics. The best-known of these equa- 
tions is Mathieu’s equation. A somewhat more complicated equation is 


d*y 
vr + [026-2 + O_1e~¥ + Oo + Ore + B2e**¥]v = 0 (1a) 
dy? 
which reduces to Mathieu's equation for 
65 = 6, 0. =0,= 61, 6» = 6 = 0, 


where the asterisk is used to denote the conjugate complex quantity. 
This paper is concerned with the determination of the solutions 


+0 
o(y) = e% DY creit¥ (2) 
of Eq. (1a) subject to the restrictions 
65 = 4, 0", = 6, 6" = 02 (1b) 
and 
0,=O(u), 62 = O(u*), (3) 


where yp is a small positive quantity. It will be seen that solution of the problem in- 
volves the determination of the “characteristic exponent” o from the equation 


sin iro = \/D sin rvV/O0, (4) 
where D denotes the expansion 
DP = 1+ C6 + Ce + Cyn + C#5? + Cade + --- (S) 


in the three real combinations 


2 2 
5 = 0_161, € = 6_fo, n = $(0:0_2 + 0_102) (6a) 


9 


of the four quantities, real and imaginary parts of 6, and 62. D is a power series in wu 
since 
§=O(u’), €=Ol(u), 1 = O(n). (6b) 


The coefficients C of the series depend on @o alone. 

The numerical evaluation of the coefficients of the expansion is the principal aim 
of this paper. This is best accomplished by first re-expressing the “doubly infinite” 
Hill determinant D in terms of its “simply infinite” principal subdeterminants D,, 


* Received June 13, 1946. 
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D = f(Do, D,, Da, - cis i, (7) 
and then expanding D, into the series 
n n n n 2 
De = 14+A8+Ac+AatAe +-:-. (8) 


The coefficients of the expansions (7) and (8) are tabulated in Tables II and I re- 
spectively for a convenient range of 69. For the sake of simpler printing the notation 


Ava = jn, Siein* } (8’) 


will be used whenever the subscript of A becomes excessively long. 

The practical solution of Eq. (1) is carried out in four steps. First, the determi- 
nants D,, Eqs. (8), are evaluated by means of Table I. Next D, Eq. (7), is determined 
from Table II. The third step consists in solving Eq. (4) or one of its variants (13a, b, 
c) for ¢, and the last step is the determination of the coefficients c, of solution (2). A 
convenient method for carrying out this last step is discussed in Section 2. The deriva- 
tions of the formulas for {m, 5‘e’n*} and for the coefficients of (7) are presented in 
Section 3. A numerical example is given in Section 4. 

The present paper is based on a study which was recently undertaken at the 
McDonnell Aircraft Corp. under the sponsorship of the Bureau of Aeronautics, U. S. 
Navy Department. The study was prompted by recent instances of control difficulties 
of some helicopters and rotor blade failures of others. As will be shown in a separate 
paper,’ the natural modes in which hinged rotor blades flap can be represented by 
solutions of Eq. (1) multiplied by suitable damping factors. It will be found that the 
stability of the blade motion decreases as the speed of advance of the helicopter in- 
creases (as uw increases). Nevertheless, instability does not set in, because an aero- 
dynamic damping effect outweighs, at all feasable speeds, the tendency towards in- 
stability which results from the flapping motion. 

The writer’s thanks are due to his colleague, Elizabeth J. Spitzer, for checking the 
derivations and the numerical work. The writer also wishes to express his indebted- 
ness to Messrs. W. R. Foote, H. Poritsky and J. J. Slade, who in their paper on 
rotational instability of shafts? applied a Laplace expansion to a doubly infinite de- 
terminant, and thus suggested the present approach. 

2. Method of solution. The solution of Eq. (1a) is assumed in the standard form 


oy) = ed) cret*¥, (2) 


Substitution of expression (2) into Eq. (1a) leads to the infinite set of homogeneous 


equations for the coefficients c,(¢): 


k= —?2: Boc- atic at [(o— 27)?+6o |c_2+6_¢_1+6_2¢0 a 0, 
k= —1: 626 s+6:c_2+ [(o—i)?+0o |c_1+6_1c0+6_2¢1 =(, (9) 
k=0: 6o¢_2+61¢_1+ [o?+-0o |co+0_161+0_2¢2 =0, 


1G. Horvay, Rotor blade flapping motion, to be published soon. 
2W. R. Foote, H. Poritsky and J. J. Slade, Critical speeds of a rotor with unequal shaft flexibilities, 


mounted in bearings of unequal flexibility, Journal of Applied Mechanics, 10, A77, 1943. 
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k=—1: B2c_1+0160+ [(o+4)?+60 |e: +0_162+0_2¢3 =0, 
k=—2: B2co+6161+ [(o+2%)?+6 }co+6_163+0_2¢4=0, 


The equations are consistent if their determinant, A(¢), vanishes. The consistency 
criterion 
A(c) = 0 (10) 
can be expressed in the much simpler form?‘ 
sin irs =. + VD sin rV/Oo, (4) 
1 612 b_2V2 0 0 | 
A1y1 1 6191 621 0 


where 


D= A(0) = . O2Vo0 410 1 0_1V0 6_2Vo (11a) 


| 
0 Ooy1 Ay 1 6191 . | 





0 0 O2V2 A1V2 1 


is the determinant of system (9) for ¢=0 when each equation is divided by the coeffi- 
cient of the diagonal term, and 


1 


a 11b 
._F (11b) 


ye = 
D is either positive or negative, and so is 0. Thus the quantity ~/D sin 7/4, is 
either real or pure imaginary. In the first case set 
g = VD sin rV/0o = — V/—D sinh rV/ — 0. (12a) 
In the second case set 
Vg =9/i= —V—Dsin rV0 = — VD sinh r\/— bo. (12b) 


Then the solution of the transcendental equation (4) is given by 





1 ae 
o = — log (q/i + V1 — q?) + mi (m = 0,+1,+2,---), 

T 
1 —4q ; 

= — arctan ——— + mi for —1Sq381, (13a) 
T V/1 —¢ 
1 chains 

= — log (g++ Vg? —1) + (m—4)i for gS —1,921, (13b) 
T 
1 a . , 

= — log (q’ + Vq’2?+ 1) + mi for g imaginary. (13c) 
T 


’ Whittaker and Watson, A course of modern analysis, Cambridge, 1927, p. 416. 
*M. J. O. Strutt, Lamésche, Mathieusche und verwandte Funktionen in Physik und Technik, Springer 


1932 (Edward Bros., 1944), p. 22. 
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Once D is known, the calculation of ¢ from (13) is simple matter. For any m there 
are two solutions a; and o2 which differ only in sign 


(14a) 


Os = — 0}. 


Let o; be the solution with the positive real part 


o1 = or + toi, o, = 0. (14b) 


The function v:(W) (or vo(W)) associated with o; (or with a2) is readily obtained by plac- 
ing a; (or g2) into the system (9) which is now limited to the equations k= —N, 
—N+1,---, —1, +1,---, +N. Assuming c,=0 for k<—WN and k>+N, one 
can solve the 2N equations for c_w, C_w41, °° - , C1, G1, * * * , €w in terms of the arbi- 
trary constant co, and then use equation k=0 asa check. The greater is N, the more 
harmonics are taken into account, and the more accurate is the solution. In practice 
the calculations are most conveniently carried out by solving equations k= N and 
k=—WN for cy and c_y in terms of the variables cy_1, cv-2 and c_n41, C_nv+42, respec- 
tively. The results are then substituted into equations k=N—1 and k= —N-+1; 
similarly cy_; and c_y4; are determined. Continuing the process one finally arrives at 
equations k= +1 and k=—1 involving the two variables ¢ and c_; only, and the 
parameter Cp which can be assumed as 1. One eliminates one of the unknowns, say c_,; 
determines from the real and imaginary parts of the remaining equation the real and 
imaginary parts of c,, and then, retracing the steps, obtains in succession the numeri- 
cal values of c_, C2, C_2, - + * , €n, C_n. 

Evidently, in principle, it is immaterial what mi (m=0, +1, +2, -- - ) is used 
in the o of Eqs. (9). For instance, the set of equations k= —N to +N with m=2, 
the set of equations k= —N—2 to +N-—2 with m=4, and the set of equations 
k=—N+3 to +N+3 with m=-—1 are identical. Thus, as one passes to the limit 
N-«, any m and any 2N+1 adjoining equations will lead to the same function 
vi(W) [or v(W)]. In practice, where one is limited to a finite number of equations, 
2N+1, it is best to use the centrally located equations k= —N to +N with an m 
which makes cp the dominating term in the series (2). 

In general 7, associated with oi, and v2, associated with o2, are linearly independent 
functions. An exceptional case arises when o,=0, and @; is an integral multiple of 4. 
Then substitution of 0; and o2 into the system (9) leads to the same function 

40 
1= ery cpe**¥, (v = 0 or $). (15a) 


=~ 


The second, linearly independent, solution is now a “quasiperiodic function” : 
+2 +2 
oy = et ly Sache 3 ine, (v = 0, }). (15b) 
—2 —< 


For convenience the functions (15a) will be called “purely periodic” functions. 
Determination of the purely periodic solutions (15a) forms the subject matter of 
most investigations on Mathieu and Hill differential equations. The purely periodic 
‘9M. J. O. Strutt, loc. cit., p. 23. As an exception there may be two purely periodic solutions. For in- 
stance, for 09=4, 6:=62.=0, one obtains 1; =cos 2y, v2=sin 2y. 
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solutions are usually of greatest interest, because they separate the y-regions of sta- 
bility (7, =0; v; and v are oscillatory) from the y-regions of instability (¢, >0, 1-7 © 
as YW). A purely periodic solution can be obtained, cf. Eqs. (12), (13), only when 
q' =0 (D=0, or perhaps 0) =k?), or when g= +1 (Dand 6¢ are such that VD sin rvV/ 85 
=1 or~/—D sinh r\/ —6,=1). In general 6, =k? does not provide a purely periodic 
function. 

In the present analysis the principal interest is attached to the unstable solutions 


+00 

(vy) = ee >> cyeieit®v, (o, > 0) (16) 

of (1) which, after multiplication by a damping factor e~*¥/?, are still stable. These 

solutions are in the “transition region” which extends from the y-value for which v() 

is purely periodic to the u-value for which e~"¥/? v() is purely periodic. It will be seen 

in Reference 1 that a rapidly advancing helicopter usually operates in the transition 
region. 

3. Expansion of Hill’s infinite determinant. It will be convenient to call the de- 
terminant D, Eq. (11a), a doubly infinite determinant to indicate that it extends to 
infinity both upward and downward. Simply infinite determinant are the principal 
subdeterminants of D, 


1 O1y¥n 9_29n 0 > | 


A1Vn41 1 O_1Yn41 O2Vn41 r cy 1 9_1Vn43 9_2Yn43 
Dn =| Oryng2 O1Yn42 «1 Oayng2 «13 Em =| * Ory «1 Laying O-19my2|. (17, b) 
0 Oxyn4a Yaga 7 ie © Oafass Crags 1 age 
- | i O bxyn 1 1 | 


The first extends to infinity downward, the second upward. Simply infinite determi- 
nants are also the auxiliary subdeterminants 


| O14" 61, O_2Vn 0 - | | . 
O29n41 1 O1y7n41 O-2¥n41 & 1 — B1Yn43 O-29n43 0 | 
S, = 0 O1Vn42 1 Sams «6+ 13: Tyme i A1Vn42 1 Oyny2 O |. (18a, b) 
0 A29n43 AVn43 1 , , O29n41  A1Vn41 1 O29n41 


le 0 O2¥n An O_1yn | 
S, differs from D, only in the first column; 7, differs from £, only in the rightmost 


column. 
One readily establishes the recurrence relations 
D, = Days — O1VnSn4t + 610_2YnVn+15n42 = EVn Vn42D nas 
+ €?ynVn41¥n42¥n4+3Dn+4s (19a) 
Sn = O1:9nDnz1 — 9-1029nVn41Dage + €VaVni1Sn-2- (19b) 


A Laplace expansion of the doubly infinite determinant D along a dividing line 
between row k=0 and k= —1 leads to the following expression involving only simply 
infinite determinants of type D,, E,, S,, Tn: 

D = DoE; — SoT1 + yoyi(6-102D1T2 + 0:0_2S1E2) 
2 
= €(yov2Di Es + yiD2E2 + yovrS1T 2) 
2 2 2 
+ € (yor yevsDi Es + YoriveE2Ds + yoviyeD2Es). (20) 
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Since, by virtue of (1b) 


E, = Dr, = Dy (21a) 
T, =S. (Sn, when 6; ¥ 01, 02 ¥ 62) (21b) 

and by virtue of (2) and (19b) 
SnTm = O(n’), (21c) 


a replacement of E, by D, and a repeated insertion of (19b) into (20) gradually elimi- 
nates all but the D, type of determinants from the expression for D. It is also found 
that 61, 0_1, 62, 6-2 appear only in the combinations 6, €, n given in (6a). Thus, by virtue 
of (6b), the expansion of D progresses in powers of u?. Using the notation 
2 
Your = Yoviye, (22a) 
one finds that to u!® terms 
2. 2 
D = DoD; — byo1D1De — €(yo2DiDs + yirDe) + n(2¥o11D2 + Z¥012D1Ds) 

— 5e(4yo112D2D3 + 2yo123D) D4) + €*(yor2sD 1D, + 2 ¥o112D2D3) 

+ en(4Vo1123D 2D + 2 Yo1234D) 1D, + 2-yo1122D3) 

—_ 5e?(4-yor1236D2D5 + 4yo11223D 3D + 2 -vo123451)1De) T-*. (23) 


The same process can also be carried out for the simply infinite determinant D,,. 
Disregarding the exceptional case 69.=k? (y,= ©), one finds that to p’® terms 


Do = Di — byoD2 + (— €Vo2 + 2nyo12)D3s + (— 25 + €*7) ores) 4 
+ 2enyorrssD5 — 25€*Yor2s45sDe +--+: , (24a) 


and D, is obtained from Dy by increasing the subscripts in the latter’s expression by n. 
(24b) 


It will be convenient to introduce at this point the notation 
p» Vase = Vase + Yasr + Vers t °°: , (22b) 
p v2 >, Yass = v2, V356 + yo3 V467 + Yeu), Tor (22c) 
Noting that 
lim D, = 1, (25) 


one obtains, by repeated application of (24), 

D=(1- 5y01)De + (- bY12 — €Yo2 + 2nyo12)Ds ele 
[1 — (yor + y12)6 — yore + 2yo12n|Ds + --- - (26a) 
1+ Ab+AetAmnt---, (26b) 


where 


wn Se 2 hes £,03F ve Oo Eimedm:--. 
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The coefficient of a general term, like 5e?, is obtained as follows: Equation (24a) 
gives rise to the following symbolic products containing de’: 


[— 5yo1][— €vo2][— €yo2], (a) 
[— dyo1] [e*yor2s], (b) 
[— eyo2|[— 2d yous], (c) 
[1] [— 26€*yorss4s]. (d) 


It is found that (a) contributes 


= Zz Yo > youd, Ja a yor >, Ysa, 7a ba yor ¥35 >, V67 (28a) 


to A$... (Note that no subscripts can be repeated, nor can any be skipped as one 
passes from one )> to the next )_; furthermore (a) gives rise to three distinct summa- 
tion expressions, because yx ,,41 Can appear in the first place, in the second place, and 
in the third place.) The relation (b) yields 


=~ y yor), Yo345 — > yous >, 455 (28b) 
2>° vor >, Ysuss + 2D, Yorss >, Yas (28c) 


— 2>0 yorssas- (28d) 


The contributions (28a, b, c, d) sum up to {0, 5e?}. By increasing the subscripts of 
the expressions (28) by 2, one obtains 


Ave =—— > ¥23 8 >, bea i yu d, Ys), 7a p yu D> ¥sr Vso 
= bi yes Do ¥4567 — > yess), Yer + > x youd, 5678 
+2 , Yosas Ves — 2); V234567- (29) 


The determination of the other {m, d‘e‘n*} is similar. 

The numerical values of the coefficients {n, dicin*} are given in Table I to 5 
decimal places, for 6. ranging from +0.9 to —1.0 (the interesting range in helicopter 
theory). In the evaluation of {n, d‘e’n*} the first 51 ym were taken into account. 
(The accuracy obtainable is thus equivalent to the use of a 101-row approximant to 
D.) yo to Yeo were computed in some instances to 6, in some instances to 7 decimal 
places; yo: to Ysoe were computed to 7 decimal places. It is expected that the entries of 
Table I are in error by not more than 2 units in the fifth decimal place.® 

It is readily seen that the present method is not limited to Eq. (1), but can be 
extended to the general Hill differential equation where 6,,e‘"¥ form a convergent 
series. 

For the special case of Mathieu's equation (e€=7=0), one finds by (23) and (27) 
that 


which (c) yields 


and (d) yields 





. 1 
D,(Do — 5yo.D2) = 1 — 28 >> 


p . 
ko 09 — k® 0) — (k + 1)? 


+ O(8?) 
sit 05 7 cot Try Oy 4 0(82) (30) 
bali » "(480 — 1)V/80 


6 An experienced computer can calculate a column of Table I in somewhat less than a day. 
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This formula was used by H. Bremekamp in 1926 in a study of the flow of electrons 


in metals. ® 

4. Example (a). Given 0=0.2, 6,;=0.19685+0.334657, 62=0.03875 —0.10258:. 
Determine 7, ve. One finds 6=0.15074, e=0.01202, » = —0.01635, and, by Table I, 
Dy =1.8422, D, =0.9434, D,=0.9934, D;=0.9980, Ds=0.999, Dj =De=1.000. Like- 
wise, by Table II, D=2.3291. Therefore, g=VD sin 1v/6,=1.5052 and by (13b) 
o,=0.30782+i/2, o2= —0.30782 —i/2. The associated functions v,(W) and w(yW) are 
determined from the equation system (9). Normalizing to co=1, and using the equa- 
tions k= —4 to k= +4, one obtains’ 


te |< 


v 3y 
v,; = (— 0.1898 + 1.9817 i)ersamee | — 0.0958 cos — + sin "1 + 0.2076 cos 4 


= 


3 5 Sy iy 
— 0.0083 sin ud — 0.0125 cos ns — 0.0038 sin — + 0.0008 cos 7 


+ 0.0019 sin = arr \ 


a 
s 


v * 3y 
(1.6652 + 0.74672) e—°-3078¥ {eos = + 0.4484 sin 9 + 0.2107 cos > 


_ sy Sy ; _ ow Ty 
+ 0.0188 sin . + 0.0041 cos = + 0.0101 sin . + 0.0005 cos ry 


~ - 


iv 
+ 0.0020 sin : = ie \ 


5. Example (b). Given 0)=0, 0; =0.37249 +0.633231, 02 = 0.13875 —0.367282. De- 
termine o. One finds 
g=\ D6o-sin ry 00/2/00 = 2v/0.4392 = 2.082, 
1} => —- 62 = 0.4339 + 1/2. 





6M. J. O. Strutt, loc. cit., p. 26. 
7 Note that the use of «;=0.30782—i/2 leads to the above expression of 1; when c; is normalized to 


1, and to the conjugate complex of the above when ¢o is normalized to 1. 
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TABLE I.* Numerical values of {n, didnt}. 









































60 = 9 8 7 -6 s | .45 4 -35 
{0, 3} | 7.83180 4.63590 3.70199 3.38325 3.38203 3.48245 3.65865 3.92977 
{0, e} : = . 90665 — .24871 -00257 - 16467 - 30899 . 38655 -47423 . 57886 
{0, »} | 6.36577 3.51924 2.63695 2.27050 2.14607 2.15140 2.20218 2 
{0, 67} — .55015 — .30117 — .22352 — .19037 —.17858 — .17822 — .18162 — .18937 
{0, de} .52473 . 27858 . 20047 . 16578 . 15052 .14789 - 14837 -15229 
{0, dn} — .15952 — .08611 — .06307 — .05310 — .04908 — .04867 — .04928 — .05105 
{0, 6} .00265 .00143 -00104 .00087 .00081 .00080 -00081 ° 
{0, e} — .41427 — .22658 —.16801 —.14320 — .13403 — .13372 — .13622 —.14199 
{0, en} | —,00590 — .00363 — .00301 — .00278 — .00280 — .00292 — .00. — .00330 
{0, 9} | —,.00638 — .00342 — .00249 — .00208 — .00191 — .00189 — .00190 — .00196 
{0, de} |} —,.00463 — .00248 — .00178 — .00150 — .00136 — .00133 — .00133 — .00136 
{0, 6} | .00049 -00025 -00019 .00015 00014 -00014 -00013 00014 
{O, de?} .00317 .00170 -00123 .00103 .00095 .00093 .00094 .00096 
{0, den} — .00021 — .00012 — .00009 — .00008 — .00007 — .00007 — .00007 — .00007 
{1, 8} | —3.27931 —1.61410 —1.05991 — .78342 — .61797 — .55795 — .50802 — .46584 
{1, e} | —1.26507 — .63934 — .43033 — .32553 — .26241 — .23943 — .22021 — .20392 
{1, 9} | —.80268 — .38700 — .24905 — .18048 — .13964 — .12488 —.11264 — .10233 
{1, 67} -04439 -02131 -01366 -00986 .00759 00678 -00610 -00553 
{1, de} | —,01646 — .00763 — .00471 — .00328 — .00243 — .00213 — .00187 — .00166 
{1, dn} .00739 -00351 00223 .00160 -00122 .00109 .00098 .00088 
{1, 6} | —.00009 — .00004 — .00003 — .00002 — .00002 — .00001 — .00001 — .00001 
{1, e} -03150 -01513 .00970 .00701 .00540 -00482 -00434 -00394 
{1, en} -00132 .00063 .00041 -00030 .00024 .00020 .00019 .00017 
{1, 77} | -00019 -00009 - 00006 .00004 .00002 .00002 .00002 .00002 
{2, 8} — .05351 — .05160 — .04981 — .04813 — .04654 — .04579 — .04505 — .04434 
{2, €} | —.03050 — .02958 — .02872 — .02791 — .02714 — .02678 — .02642 — .02607 
{2, 9} | —.00619 — .00591 — .00565 — .00541 — .00519 — .00508 — .00498 — .00487 
{2, 8} | -00025 -00024 .00023 -00021 -00021 .00020 -00020 .00019 
{2, de} -00001 -00001 -00001 -00001 -00001 -00001 -00002 .00002 
{2, nm} -00002 -00002 -00002 -00002 .00002 .00002 .00002 -00002 
{2, €} .00017 -00017 -00016 .00015 .00015 .00014 .00014 -00013 
{3, 3} | —.01368 — .01349 — .01330 — .01311 — .01293 — .01284 — .01275 — .01267 
{3, €} | —.00914 — .00903 — .00892 — .00881 — .00871 - 866 — .00861 — .00856 
{3,7 |} —,00092 — .00090 — .00088 — .00086 — .00086 — .00084 — .00084 — .00083 
{3, 6} 00003 00003 00003 .00003 .00003 00003 00003 00003 
{3, de} | 00001 00001 00001 .00001 .00001 00001 00001 00001 
{3, e] | 00002 0000 2 0000 2 .00002 -00002 00002 00002 00002 
{4, 6} — .00551 — .00546 — .00543 — .00538 — .00534 — .00532 — .00530 — .00528 
{4, €} — .00401 — .00399 — .00396 — .00393 — .00391 — .00390 — .00388 — .00387 
{4, n} — .00024 — .00024 — .00023 — .00023 — .00023 — .00023 — .00023 — .00023 
{4, de} .00001 00001 00001 00001 .00001 00001 00001 00001 
{5, 5} — .00276 — .00275 — .00273 — .00272 — .00271 — .00270 — .00269 — .00269 
{5, e} — .00213 — .00212 — .00211 — .00210 — .00209 — .00209 — .00208 — .00208 
{5. 9} | —,.00008 — .00008 — .00008 — .00008 — .00008 00008 — .00008 -- .00008 
{6, 6] | —,00158 — .00157 — .00157 — .00156 — .00155 — .00155 — .00155 — .00155 
{6, €} — .00126 — .00126 — .00126 — .00125 — .00125 — .00125 — .00125 — .00125 
{6, 9} — .00003 — .00003 — .00003 — .00003 — .00003 — .00003 — .00003 — .00003 
Oe = 3 -25 o” -15 1 -05 o** 
{0, 5} 4.33216 4.93480 5.87874 7.49588 10.78515 20 .74569 
{0, e} - 71097 - 88889 1.14867 1.57391 2.41481 4.92154 . 250000 
{0, »} 2.48046 2.75850 3.21013 4.00081 5.62957 10 .59568 . 
{0, 6} | —.20279 — .22456 — .26021 — .32296 — .45255 — .84825 — .039865 
{0, de} - 16054 -17498 . 19957 . 24378 - 33621 -62017 -028683 
{0, dn} | —,.05432 — .05977 — .06883 — .08490 —.11824 — .22028 — .010290 
{0, 5} | .00089 -00098 -00112 -00138 -00192 -00357 -000166 
{0, e} |} —.15201 — .16827 — .19494 — .24186 — .33884 — .63496 — .029835 
{0, en} | — .00364 — .00416 — .00496 — .00630 — .00905 — .01736 — .000836 
{0, »*} |; —,.00208 — .00228 — .00262 — .00322 — .00448 — 00832 — .000388 
{0, de} — .00144 — .00158 — .00180 — .00221 — .00305 — .00566 — .000262 
{0, 5} | .00014 -00015 -00018 .00022 -00032 .00059 .00002 
{0, de*} } 00103 -00112 .00129 -00160 -00222 -00412 -000194 
{0, den} | —.00007 — .00007 — .00008 — .00009 — .00012 — .00023 — .000012 
{1, 3} — .42975 — .39853 — .37126 — .34726 — .32596 — .30694 — .28987 
{1, e} — .18993 —.17778 — .16712 — .15769 —.14929 — .14176 — .13496 
{1, 9} — .09354 — .08595 — .07934 — .07355 — .06844 — .06388 — .05980 
{1, 8} .00505 -00463 -00426 -00395 -00366 -00341 -00319 
{1, de} | —,00149 — .00134 — .00120 — .00109 — .00099 — .00091 — .00083 
{1, én} -00080 -00073 -00067 -00062 -00057 -00053 -00050 
{1, 8} — .00001 — .00001 — .00001 — .00001 — .00001 — .00001 — .00001 
{1, €} -00359 -00330 -00303 -00281 -00261 -00243 -00227 
{1, en} -00014 -00014 -00012 -00012 -00011 -00011 -00010 
{1, 9} -00002 -00002 -00001 -00001 -00001 -00001 -00001 








* Note 1. {, d%e&n*} which are less than .00001 in magnitude, or for which » >6, are not shown. 
0 0 
& Ae De 
** Note 2. For @o=0 (yo= @) the coefficients * .*** of — are given. 
ye ye ye 
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TABLE [. (Continued) 

6 = | 25 a 15 a 05 ore 

{2, 5} — .04365 — .04297 — .04232 — .04168 — .04106 — .04045 —.03987 
{2, e} — .02573 — .02540 — .02507 — .02475 — .02445 — .02415 — .02385 
{2, 9} — .00477 — .00468 — .00458 — .00450 — .00442 — .00433 — .00425 
{2, 8%} .00019 .00018 .00018 .00018 .00017 .00017 .00017 
{2, de} -00002 .00002 .00002 .00002 .00002 .00002 -00002 
{2, dn} .00002 .00002 .00002 .00002 .00002 .00002 .00002 
{2, e} -00013 -00013 -00013 .00013 .00012 00012 .00012 
{3, 8} — .01258 — .01250 —.01241 — .01233 — .01225 —.01217 — .01209 
{3, €} — .00851 — .00846 — .00842 — .00837 — .00832 — .00827 — .00823 
{3, n} — .00082 — .00081 — .00080 — .00079 — .00078 — .00078 — .00077 
{3, 87} .00003 .00002 .00002 .00002 .00002 .00002 .00002 
{3, de} .00001 .00001 .00001 .00001 .00001 .00001 .00001 
{3, e2} -00002 00002 .00002 .00002 .00002 .00002 .00002 
{4, 5} — .00526 — .00524 — .00522 — .00520 — .00518 — .00516 — .00514 
{4, €} — .00386 — .00385 — .00383 — .00382 — .00381 — .00380 — .00378 
{4, »} — .00022 — .00022 — .00022 — .00022 — .00022 — .00022 — .00022 
{5, 5} — .00268 — .00268 — .00267 — .00266 — .00266 — .00265 — .00264 
{5, ¢} — .00207 — .00207 — .00207 — .00206 — .00206 — .00205 — .00205 
{5, »} — .00008 — .00008 — .00008 — .00008 — .00008 — .00008 — .00008 
{6, 5} — .00155 — .00155 — .00154 — 00154 — .00154 — .00153 — .00153 
{6, €} —.00124 — .00124 — 00124 — .00024 — .00124 — .00123 — .00123 
{6, n} — .00003 — .00003 — .00003 — .00003 — .00003 — .00003 — .00003 

0 = —.05 —.1 —.15 —.2 —.25 —.3 —.35 —.4 

{0, 5} —19.32207 —9,35137 —6.04482 —4.40273 —3.42538 —2.77963 —2.32283 —1.98371 
{0, e} —5.06707 —2.56220 —1.72447 —1.30380 —1.05015 — .88013 — .75802 — .66590 
{0, »} —9.46238 —4.48740 —2.84360 -—2.03119 —1.55045 —1.23488  —1.01322 — .84989 
{0, 6} .75141 .35491 . 22402 15940 .12120 .09616 .07860 .06569 
{0, de} — .53194 —.24724 —.15352 —.10745 — .08037 — .06274 — .05042 — .04144 
{0, én} .19281 09051 .05679 .04017 .03036 .02396 .01947 .01617 
{0, 63} —.00311 ~—.00146 — .00092 — .00065 — .00048 — .00038 — .00031 — .00026 
{0, e} -56222 i .26550 -16755 -11919 .09061 .07188 .05874 .04908 
{0, en} .01609 .00775 .00500 .00363 .00281 .00226 .00188 .00160 
{0, n?} -00724 .00338 .00211 .00149 -00112 .00088 .00072 .00059 
{0, 5%} .00488 .00227 .00141 .00098 .00074 .00059 .00047 .00039 
{0, 8} ~ .00048 — .00023 — .00014 — .00010 — .00008 — .00006 — .00006 — .00004 
{0, de?} — .00357 — .00167 — .00103 — .00073 — .00055 — .00044 — .00035 — .00029 
{0. den} .00023 .00010 .00006 .00004 .00003 .00003 .00002 .00001 
{1, 5} — .27445 — .26046 —.24772 — .23607 — .22538 — .21553 — .20643 — .19800 
{1, e} — .12880 —.12318 —.11805 — .11332 — .10897 —.10494 —.10120 — .09772 
{1,9} — .05614 — .05282 — .04981 — .04707 — .04457 — .04227 — .04016 — .03821 
{1, 67} .00299 .00281 .00264 .00249 .00235 00223 .00211 .00201 
{1, de} — .00076 — .00069 — .00064 — .00059 — .00054 — .00050 — .00047 — .00043 
{1, 5} .00046 .00043 .00040 .00038 .00036 .00034 .00032 .00030 
{1, 3} — .00001 — .00001 — .00001 — .00001 

{1, e} .00213 .00200 .00188 .00178 .00168 .00159 .00151 .00143 
{1, en} .00009 .00009 .00008 .00008 .00007 .00007 .00007 .00006 
{1, 97} 00001 .00001 .00001 .00001 .00001 .00001 .00001 

{2, 5} — .03929 — .03873 — .03819 — .03766 — .03714 — .03663 — .03614 — .03566 
{2, e} — .02357 — .02328 — .02301 — .02275 — .02249 — .02223 — .02198 — .02174 
{2, 7} — .00417 — .00409 — .00401 — .00394 — .00387 — .00380 — .00373 — .00367 
{2, 87} .00016 .00016 .00016 .00015 .00015 .00015 .00015 .00014 
{2, de} .00002 .00002 .00002 .00002 .00002 .00002 .00002 .00002 
(2, 6n .00002 .00002 .00001 .00001 .00091 .00001 .00001 .00001 
{2, e} .00011 .00011 -00011 .00011 .00011 .00010 .00010 .00010 
{3, 5} —.01201 — .01193 — .01185 — .01178 — .01170 — .01163 —.01155 — .01148 
(3, e} — .00818 — .00814 — .00809 — .00805 — .00801 — .00796 — .00792 — .00788 
{3, 7} — .00077 — .00076 — .00075 — .00075 — .00074 — .00073 — .00073 — .00072 
{3, 87} .00002 .00002 .00002 .00002 .00002 .00002 .00002 .00002 
{3, de} .00001 .00001 .00001 .00001 .00001 .00001 .00001 .00001 
{3, e} .00002 .00002 -00001 .00001 .00001 .00001 .00001 .00001 
{4, 3} — .00512 — .00510 — .00509 — .00507 — .00505 — .00503 — .00501 — .00500 
{4, €} — .00377 — .00376 — .00375 — .00374 — .00372 — .00371 — .00370 — .00369 
{4, 7} — .00022 — .00022 — .00021 — .00021 — .00021 — .00021 — .00021 — .00021 
{5, 3} — .00264 — .00263 — .00262 — .00262 — .00261 — .00261 — .00260 — .00259 
{5, €} — .00204 — .00204 — .00204 — .00203 — .00203 — .00202 — .00202 — .00202 
{5, n} — .00008 — .00008 — .00008 — .00008 — .00008 — .00008 — .00008 — .00008 
{6, 3} — .00153 — .00153 — .00152 — .00152 — .00152 — .00152 — .00151 —.00151 
{6, « — .00123 — .00123 — .00123 — .00122 — .00122 — .00122 — .00122 — .00122 
{6, 7} — .00003 — .00003 — .00003 — .00003 — .00003 — .00003 — .00003 — .00003 
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TaBLE I. (Continued) 
—.45 —.5 —.55 —.6 -.7 -.8 -.9 —1.0 
—1.72274 —1.51622 1.34910 1.21138 — .99849 — .84237 —.72362 — .63067 
— 59385 — .53588 — .48819 — .44823 — .38497 — .33706 —.29947 —.26917 
—.72520 — .62733 —.54879 — .44863 — .38670 —.31617 — .26349 —.22300 
05584 04813 04195 .03691 02924 02374 01965 01652 
— .03464 — .02936 — .02516 — .02176 — .01666 — .01306 — .01043 — .00846 
01368 01172 .01016 00888 .00695 00559 00457 00380 
— .00022 —.00019 — .00016 — .00014 —.00011 — .00009 — .00007 - 
04172 03595 .03133 02757 02184 01772 01467 .01233 
00139 00121 00107 .00096 00078 00064 00055 .00047 
.00050 00043 00037 00032 00025 .00020 .00016 00013 
00033 00028 00024 00021 .00016 00013 00010 .00008 
— 00003 — .00003 — .00002 — .00002 — .00001 — .00001 — .00001 — .00001 
— .00025 —.00021 — .00018 — 00016 — .00012 — .00010 — .00008 — .00007 
00001 .00001 .00001 .00001 
—.19017 — 18288 —.17608 —.16971 —.15815 —.14793 — .13882 —.13067 
—.09448 —.09144 — .08859 — .08591 — .08102 — .07664 — .07272 — .06917 
—.03641 — .03473 — .03318 — .03173 —.02911 — .02682 — .02480 — .02300 
00191 00182 00173 .00166 00151 00139 00128 00118 
—.00040 — .00037 — 00034 — .00032 — .00028 — .00024 — .00020 — .00018 
:00029 .00027 .00026 00025 00022 .00020 00018 00017 
.00136 00130 00124 00118 .00108 .00099 .00092 .00085 
00006 .00006 00006 .00006 00005 .00005 -00004 .00003 
— 03519 — .03473 — .03428 — .03384 — .03300 — .03219 —.03141 — .03067 
02150 —.02126 — .02103 —.02081 — .02037 — .01995 — .01955 —.01917 
— .00361 —.00354 — .00349 — 00342 — 00331 — .00320 — .00310 — .00150 
00014 00014 00014 00013 00013 00012 00012 00012 
00002 00002 00002 .00002 00002 .00002 .00002 00002 
:00001 00001 .00001 .00001 .00001 .00001 .00001 .00001 
:00010 00010 00010 00009 .00009 .00009 .00008 .00008 
—.01141 —.01134 —.01127 — .01120 — .01106 — .01093 — .01080 — .01067 
— .00784 —.00779 — .00775 —.00771 — .00763 — .00755 — .00748 — .00740 
— 00072 —.00071 —.00071 — 00070 — .00068 — .00067 — 00066 — .00065 
00002 00002 00002 00002 00002 .00002 .00002 .00002 
00001 00001 .00001 00001 .00001 00001 .00001 .00001 
00001 -00001 .00001 .00001 00001 00001 .00001 .00001 
— 00498 — 00496 — 00494 — .00492 — .00489 — .00485 — .00482 — .00479 
— 00368 — .00366 — 00365 — .00364 — .00362 — .00360 — .00358 — .00356 
—.00021 — .00020 — .00020 — .00020 — .00020 — .00020 — .00020 — .00020 
—.00259 —.00259 — .00258 —.00257 — .00256 — 00255 — .00253 — .00252 
—.00201 — 00201 — .00200 — .00200 — .00199 — .00198 — .00197 — .00197 
— 00008 — 00008 — .00008 —.00008 — .00008 — .00007 — .00007 — .00007 
—.00151 —.00151 —.00150 — .00150 — 00150 —.00149 — .00149 — .00149 
—:00122 —.00121 —.00121 —.00121 —.00121 — .00120 —.00120 — .00120 
— .00003 — .00003 — .00003 — .00003 — .00003 - 3 — .00003 
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TABLE II. Expansion of D (row 1 Xrow 2 Xrow 3, see Eq. 23). 
Numerical tabulation of the coefficients in row 3. 














row 1 1 6 € ” be 

Oo row 2 { DoD; D:D: } DiDs | Di DD; DDs | 

row 3 —Yyor —Yor —yu 2 you 2 te —4Yo12 —2 Yoiss 
oS 1 4.00000 -57143 —4.00000 16.00000 2.28571 9.14285 . 26891 
-45 i 4.04040 -62598 -—3.30579 14.69238 2.27628 8.27740 - 26623 
4 i 4.16667 -69445 -—2.77778 13.88889 2.31482 7.71606 - 26916 
.35 i 4.39561 .78278 —2.36687 13.52495 2.40856 7.41094 . 27845 
a 1 4.76190 -90090 —2.04082 13.60543 2.57400 7.35428 . 29586 
.25 | 1 5.33333 1.06667 —1.77778  14.22222 2.84445 7.58519 .32508 
Py | 1 6.25000 1.31579 —1.56250 15.62500 3.28948 8.22369 . 37380 
Ay 1 7.84314 1.73160 —1.38408 18.45445 4.07436 9.58673 -46037 
1 1 11.11111 2.56410 —1.23457 24.69135 5.69800 12.66222 -64023 
.05 1 21.05263 §.06329 —1.10803 44.32135 10.65956 22.44119 1.19100 
o* 0 1.00000 . 25000 0 2.00000 . 50000 1. 05556 
—.05 1 —19.04762 --4.93827 — .90703 —36.28118 -9.40624 —17.91666 —1.03936 
—.10 | 1 —9.09091 —2.43902 — .82645 —16.52891 -—4.43458 -—8.06287 — .48732 
—.15 1 —5.79710 —1.60643 — .75644 —10.08191 -—2.79379 -—4.85876 — .30533 
—.2 1 —4.16667 -—1.19048 — .69444 —6.94444 -—1.98412 -—3.30686 — .21566 
—.25 | 1 —3.20000 — .94118 — .64000 -—5.12000 -1.50588 -2.40941 — .16280 
—.3 | i —2.56410 — .77519 — .59171 —3.94477 -—1.19260 -—1.83478 — .12823 
— .35 1 —2.11640 — .65681 — .54870 —3.13541 — .97306 -—1.44157 — .10407 
—.4 | i —1.78571 — .56818 — .51020 -—2.55102 — .81169 -—1.15955 — .08635 
—.45 | 1 —1.53256 — .49938 — .47562 -—2.11388 — .68880 — .95006 — .07289 
—.5 1 —1.33333 — .44444 — .44444 -—1.77778 — .59259 — .79012 — .06238 
—.55 1 —1.17302 — .39960 — .41623 -—1.51357 — .51561 — .66530 — .05399 
—.6 1 —1.04167 — .36232 — .39063 -—1.30208 — .45290 — .56612 — .04718 
-.7 | 1 — .84034 — .30395 — .34602 — .98863 — .35759 — .42069 — .03687 
—.8 | 1 — .69444 — .26042 — .30864 — .77160 — .28935 — .32150 — .02953 
-.9 1 — .58480 — .22676 — .27701 — .61557 — .23869 — .25126 — .02411 
—1.0 | 1 — .50000 — .20000 — .25000 -—. — .20000 — .20000 — .02000 

row 1 é en be? 

6 row 2 bata pro HF pean D:Ds - D:D. Did. am 

row 3 Yoies 2 Yas 4 Yours 2 Yorrs 2 Yorrss —4 Yours = —4 Youres —2 Josss 

me | — .13445 —4.57143 1.07563 -01735 1.30612 -06941 -30732 .00071 
.45 — 13312 — 4.13870 96812 -01712 1.16583 -06226 .27271 .00070 
-4 — .13458 —3.85803 -89721 .01726 1.07168 -05752 - 24922 .00070 
.35 — .13922 —3.70547 - 85676 -01779 1.01520 -05475 - 23498 -00072 
on | — .14793 —3.67714 - 84532 -01885 -99382 -05384 . 22846 .00076 
oan | — .16254 —3.79260 . 86688 -02064 1.01136 -05504 -23117 .00083 
on | — . 18690 —4.11185 -93450 -02366 1.08206 -05915 - 24592 -00095 
15 | — .23019 —4.79337 1.08323 -02904 1.24503 -06834 - 28136 -00117 
ol —.32011 —6.33111 1.42273 -04027 1.62336 -08948 - 36480 -00162 
.05 —.59550 -—11.22059 2.50737 -07468 2.84066 -15722 -63478 .00300 
o* — .02778 — .50000 -11111 -00347 . 12500 -00694 -02778 .00014 
-.0S | 51968 8.95833 —1.97974 — 06476 —2.21193 — .12335 — .48882 — .00260 
—.1 24366 4.03144 — .88604 — .03026 — .98328 — .05503 — .21610 — .00121 
—.15 15267 2.42938 — .53101 — .01891 — .58539 — .03288 — .12796 — .00075 
—.2 10783 1.65343 — .35943 — .01331 — .39367 — .02218 — .08558 — .00053 
—.25 -08140 1.20470 — .26048 — .01002 — .28346 — .01603 — .06129 — .00040 
—.3 .06412 -91739 —.19728 — .00787 — .21335 — .01210 — .04588 — .00031 
—.35 05203 . 72078 —.15418 — .00637 —.16570 — .00943 — .03544 — .00025 
—.4 04318 .57978 — .12336 — .00527 — .13177 — .00752 — .02804 — .00021 
—.45 03644 .47503 —.10054 — .00443 — .10675 — .00611 — .02259 — .00017 
—.5 03119 . 39506 — .08317 — .00378 — .08779 — .00504 — .01848 — .00015 
—.55 02700 -33265 — .06967 — .00326 —.07311 — .00421 — .01531 — .00013 
—.6 02359 283 — .05897 — .00284 — .06154 — .00355 — .01282 — .00011 
—.7 -01843 21035 — .04337 — .00221 — .04475 — .00260 — .00923 — .00009 
—.8 -01476 16075 — .03281 — .00176 — .03349 — .00195 — .00684 — .00007 
-.9 -01206 12563 — .02538 — .00143 — .02564 — .00150 — .00518 — .00006 
—1.0 .01000 - 10000 — .02000 — .00118 — .02000 — .00118 — .00400 — .00005 





* Note 1. For 6:=0 (yo= «) the coefficients of Di. are given. 
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ON THE MECHANICAL BEHAVIOUR OF METALS IN THE 
STRAIN-HARDENING RANGE* 


BY 


G. H. HANDELMAN, C. C. LIN anp W. PRAGER 
Brown University 


1. Introduction. The present paper is concerned with certain stress-strain rela- 
tions purporting to describe the mechanical behaviour of quasi-isotropic metals in 
the strain-hardening range. As a preparation for a more precise characterization of 
these relations, let us consider the tension test of a metal like copper or aluminum 
which does not flow under a constant stress, but exhibits strain hardening. If the test 
involves loading only, i.e., if the reduced tensile stress' o or the tensile strain € in- 
crease throughout the test, the resulting diagram of reduced stress versus strain will 
have the general appearance of the curve OPQ in Fig. 1. On the other hand, if the 
test specimen is unloaded after a cer- 
tain point, such as P, has been reached co Q 
aléng this curve, the stress-strain dia- : 
gram for unloading is found to be very 
nearly a straight line PA which is par- 
allel to the tangent of the curve OPQ 
at O. After complete unloading, the ° 
specimen shows a permanent extension 
which corresponds to the permanent 
strain represented by OA. 

To simplify the discussion, let us p' 
assume at present that the material is 
incompressible. A longitudinal extension 
e of the isotropic specimen is then ac- 
companied by a uniform lateral con- 
traction of the magnitude ¢€/2. If the discussion is restricted to states of stress and 
strain which can be reached by a single loading followed by one complete or partial 
unloading at the most, the mechanical behaviour of the material in simple tension is 
therefore completely defined by the curve OPQ. It will be assumed in the following 
that for the materials under consideration the stress-strain diagram in simple com- 
pression (OP’Q’ in Fig. 1) is obtained by reflecting the curve OPQ with respect to 
the origin O, and that the practically important portion of the curve Q’OQ, i.e., the 
portion corresponding to small and moderate strains, is represented with sufficient 
accuracy by a development of the form 


e=o0+a30*+ao0°+:--, (1) 


where ds, @s, - - - are constants. (The coefficient of the linear term on the right-hand 
side of (1) must be unity since ¢ is the reduced stress. No even powers of o can occur 


See 








Fic. 1. Typical curve of reduced stress vs. strain. 





* Received September 17, 1946. 
1 The reduced stress is defined as the quotient of the stress by Young's modulus, 
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on the right-hand side of (1), because the stress-strain diagrams for tension and com- 


pression are assumed to be congruent.) 
In the case of simple tension or compression, the mechanical behaviour of the ma- 


terial during the first loading is readily represented by a finite relation of the form (1); 
the behaviour during the first unloading, however, is most naturally represented by 
the differential stress-strain relation 

de = da, (2) 


for this form avoids explicit reference to the state of stress at which the unloading 
began. Accordingly, it is often convenient to write Eq. (1), too, in differential form: 


de = a(o)de. (3) 
Here, a(c) =de/do = 1+3a3;07+5a;0'+ - - - equals the quotient of Young’s modulus 
by the so-called tangent modulus. To arrive at a complete analytical description of the 


mechanical behaviour of the material in simple tension and compression, we must sup- 
plement the preceding equations by analytical criteria for loading and unloading. 
For tension (¢ >0) loading corresponds to do >0 and unloading to do <0; for compres- 
sion (¢ <0) these criteria must be reversed. A satisfactory criterion for loading and 
unloading is therefore furnished by the sign of oda =d(}o"). 

The present paper is concerned with the extension of this analysis to general 
states of stress and strain which can be reached by a single loading followed at most 
by one complete or partial unloading. In the case of simple tension or compression, a 
differential stress-strain relation of the form (3) which is valid for the first loading can 
always be integrated under the initial condition e=0 for o =0 and is thus equivalent 
to a finite stress-strain relation. For more general states of stress, however, a suitably 
generalized form of the differential stress-strain relation (3) may be integrable or not. 
The distinction between differential and finite stress-strain relations for the first load- 
ing is therefore no longer a purely formal matter, but acquires physical significance. 
One of the main results of the following discussion consists in the remark that the 
assumption of a finite stress-strain relation for the first loading is incompattble with 
certain postulates concerning the mechanical behaviour under those changes of stress 
which constitute neither loading nor unloading. This is shown in Section 3. Sections 2 
and 4 are devoted to the discussion of finite and differential stress-strain relations, 
respectively. Section 5 gives a method of correlating experimental results with the 
present theory. Finally, Section 6 contains a discussion of the limitations of the 
theory. 

2. Finite stress-strain relations. Using rectangular Cartesian coordinates xj, 
(t=1, 2, 3), we denote the displacement from the standard state by u,, the strain by 
€;; and the reduced stress by o;;. For the small deformations to which the following 
discussion is restricted, the strain €;; is given by 

€ij = 3(Ui,5 + 4;,:), (4) 


where u;,; stands for du;/0x;, etc. Adopting the usual summation convention regard- 
ing repeated subscripts, we define the mean normal strain as 
¢é= Fé, (5) 


and the strain deviation as 
Cig = €4j — 0055, (6) 
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where 6,; is the Kronecker delta. Similarly, the reduced mean normal stress s and the 
deviation s;; of the reduced stress are defined as 


S= jou (7) 


and 
Sig = OG — $633. e (8) 


According to the definitions of the deviations e;; and s;;, we have 
ex = 0, Six = 0. (9) 


The task of generalizing the finite stress-strain relation (1) is simplified by the 
remark that the first term on the right-hand side represents that part of the total 
strain € which is recovered upon complete unloading. The remaining terms on the right- 
hand side of (1) accordingly represent the permanent strain. In Fig. 1 the total strain 
is represented by the segment OB, the recoverable strain by AB, and the permanent 
strain by OA. 

Setting 


| a 6 j + ijn (10) 


whére «€, denotes the recoverable and «jj the permanent strain, we may assume that 
the recoverable strain is related to the reduced stress by means of the generalized law 
of Hooke: 


é:5 = (1 + v)si; + (1 — 2n)58;;. (11) 


Here v denotes Poisson's ratio. We are then left with the task of supplementing (11) 
by a relation which expresses the permanent strain occurring during the first loading 
in terms of the reduced stress. For an isotropic material, this relation can only con- 
tain scalar constants in addition to the tensors ¢€jj, 0;; and 6,;;, and their invariants. 
Furthermore, the principal axes of ¢{/ and ¢;; must coincide. Under the pressures com- 
monly encountered in the testing of materials, no permanent change of volume is 
observed, i.e., €/ =0 and ¢/=e//. A state of hydrostatic pressure therefore does not 
produce any permanent strain, and two states of stress which differ only by a state of 
hydrostatic pressure may be expected to produce identical permanent strains. The 
permanent strain ¢€/ is thus independent of s and depends only on the deviation s;;. 
Furthermore, if the stress-strain diagrams for simple tension and simple compression 
are congruent, a reversal of the signs of all stresses may be expected to produce a mere 
reversal of the signs of all principal strains. Finally, if the ratios of the principal 
stresses are kept constant during the loading process, the ratios of the principal per- 
manent strains, too, can be expected to remain constant. 

In a recent paper,? W. Prager established the most general stress-strain relation 
which is compatible with the preceding postulates. With the notations 


J, = 35:58 ji, Ts = 4S: j8 juSkiy (12) 


and 
tig = SieSez — FJ 0bi:5, (13) 


Prager’s stress-strain relation can be written in the form 


2,W. Prager, Strain-hardening under combined stresses, }. Appl. Phys. 16, 837-840 (1945). 
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di, = Ps, TIPU s, Falta + Os, Tao tas), (14) 


where P and Q must be homogeneous in the components of the stress deviation, the 
degree of P exceeding that of Q by 4. The expressions (12) are second and third order 
invariants of the stress deviation s;; (the first order invariant s;; vanishes). The tensor 
(13) is the deviation of the square s,s,; of the stress deviation s;;. 

Combining (11) and (14), we obtain the desired generalization of the finite stress- 
strain relation (1): 


é;= (1 + v)sij + (1 _ 2v)s6;; + F(J2, Js) [PUs, Somes + Os, Iss ti]. (15) 


3. Neutral changes of stress. Inadmissibility of finite stress-strain relations. In 
the case of simple tension or compression the sign of eda =d(}0*) proved to be a satis- 
factory criterion for loading and unloading. Accordingly, one might consider the pos- 
sibility of using the sign of o;do;; as a criterion in the general case. If, however, the 
term “loading” is reserved for such changes of stress which are accompanied by a 
change of the permanent strain, this criterion is not satisfactory. Indeed, on account 
of (8) and the second Eq. (9), we have 


oi;doi; = (Si; + 56;;)(dsi; + dsé;;) = $i jdsi; + 3sds. (16) 


If loading were to correspond to o;;do;;>0, a change of stress for which ds;;=0 might 
therefore constitute loading in spite of the fact that such a change of stress is not 
accompanied by a change of the permanent strain. To avoid this difficulty, we shall 
use the sign of s;;ds;;=dJ2 as the desired criterion, an increase of J; corresponding to 
loading, a decrease to unloading. 

Whereas for uniaxial stress any change of stress constitutes either loading or un- 
loading, we have three kinds of change of stress in the general case, according to 
whether J2 increases, remains constant, or decreases. An infinitesimal change of stress 
for which dJ,=0, will be called a neutral change of stress. For instance, any change of 
stress which affects only the mean normal stress, but leaves the stress deviation un- 
touched, is a neutral change of stress. A more interesting example of a neutral change 
of stress is given by 


(co 0 0) (0 dr 0) 
| | 

vij= [0 0 0}, doi;= |dr 0 0}. (17) 
10 0 0} 10 0 0) 


Equation (17) represents the stress system which arises from a combined tension and 
torsion test of a thin walled circular cylinder. Specifically, consider such a test piece 
which is pulled to an arbitrary tensile stress o. If the traction is then kept constant 
and a small torque applied, the resulting systems of stress and increments of stress 
are represented by Eq. (17). 

Let us now suppose that for the first loading (dJ2>0) we have the finite stress- 
strain relation (15) and for unloading (dJ:<0) the generalized-law of Hooke in the 
differential form 

dé;; => (1 + v)ds;; + (1 = 2v)ds6; ;. (18) 
The simultaneous use of the stress-strain relations (15) and (18) will lead to obvious 
difficulties, unless these relations give identical strain increments for neutral changes 
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of stress. We shall show that, in general, this continuity condition is not fulfilled. In- 
deed, if (15) is written in differential form, the first two terms on the right-hand side 
equal the right-hand side of (18); the continuity condition therefore requires the van- 
ishing of the remaining terms on the right-hand side of the differential form of Eq. 
(15). Consider, for example, the stress and increment of stress given by Eq. (17). A 
simple computation will show that 











2 0 O zo? 0 0 | 
‘3; = | 0 —to 0 ly t;= |0 -—4o O I 
10 8 te) 0 0 —get} 
9 dr 0) (0 todr 0) 
ds;; = (dr 0 0| dt;; = |4odr 0 0| 
10 oO 0) 0 0 0) 


as well as dJz=dJ;=0. For this special case the differential form of Eq. (14) reduces to 
de,’ = F(J2, Js) [P(J2, Ja)dsi; + OSs, Js)Js dtis), (19) 


where Je and J; are evaluated for an arbitrary state of pure tension. Since this state 
of stress satisfies the condition for a neutral change of stress (dJ,=0), de{ must vanish. 
We find then, upon substituting the values of ds;; and dt;; previously computed, that 


F(Js, Js) [P(UJ2, Js) + 402, Ja)Js o] = 0. (20) 


Now let us return to the finite stress-strain relation, Eq. (15), for the case of pure 
tension. The first component of the strain tensor (the other non-vanishing terms 
differ from this only by a constant factor) becomes 


en = o + 3F(Js, Js) [PUs, Jao + 402, Ja)Jse |. (21) 


The invariants appearing in Eq. (20) have been evaluated for an arbitrary state of 
pure tension. Consequently, Eq. (20) is valid for pure tension and the second term in 
Eq. (21) equals zero. (A similar remark holds true for each of the other non-vanishing 
strain components.) Therefore, the stress-strain relation will reduce to Hooke’s law 
for pure tension if the continuity condition is to be fulfilled. On the other hand, we 
have seen in Section 1 that the stress-strain law for tension need not be linear. Thus 
the most general finite stress-strain law coupled with Hooke’s law for unloading will 
not be sufficiently flexible to represent a tensile test if the continuity condition is to 
be fulfilled. It is necessary, therefore, to turn to differential stress-strain relations if 
both loading and unloading are to be adequately represented. 

4. Differential stress-strain relations. A system of differential stress-strain rela- 
tions can be obtained from the properties discussed in Section 2 provided certain of 
these are rewritten in such a way as to be directly applicable in differential form. 
We shall assume that given the components of the stress tensor ¢;; and the increments 
de;; there correspond unique strain increments de,;;. This implies that the increment 
in strain, de,;;, depends only on the state of stress at the given instant, o,;, and the 
increment in stress, do,;, and is independent of the way in which this state of stress 
has been achieved provided only loading has taken place. In particular, we shall 
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assume that this dependence is such that the increments in strain are linear functions 
of the increments in stress. Thus de;; can be written in the form 


dei; = (1 + v)dsi; + (1 —- 2v)ds6;; + Cijkl dex, (22) 


where the fourth order tensor ¢;j,; is a function of o;;only. For unloading, the material 
is assumed to satisfy the differential form of Hooke’s law given in Eq. (18). Loading is 
supposed to take place when dJz.>0 and unloading occurs for d/;<0. For a neutral 
change of stress, dJ2=0, the continuity condition requires that Eqs. (18) and (22) 
coincide. Consequently, 
Cijxt don. = O whenever dJ2 = sx, ds, = 0. 

Since s,; is a deviator, dJz may also be written in the form dJz=s,; dax:. Thus the 
linear form in doy), €:;¢2 dox%1, must vanish whenever s;,; do, vanishes. The coefficients 
of do; in the two forms must be proportional or 


Cijkl = CisSar, 
where the second order tensor C;; is a function of o,; alone. The stress-strain relations 


then become 


de;; = (1 + v)ds;; + (1 — 2v)dsé;; + CijdJ2, when dJz2 0;) (23 
f is 


(1 + v)ds;; + (1 — 2v)ds6;;, when dJ2 S 0. 


dé;; = 


In a certain sense, the term C;; measures the permanent deformation. Indeed, let 
us consider the infinitesimal cycle of stress which results when first do;; is applied 
and then —da;;. We assume, in addition, that the material is being loaded when do;; 
is applied. The permanent increment in strain dé will then be 


dey; = Ci; dJ>. (24) 


Since the permanent strain is independent of a state of hydrostatic stress for pressures 
within the range normally encountered in testing of materials, the tensor C;; can only 
be a function of the components of the stress deviator rather than the stress tensor 
itself. Furthermore, there can be no permanent change in volume; that is, de; =0 or 
C;;=0. Since the tensors de;;, ds;;, and 6;; are symmetric, C;; will also be symmetric. 
In addition, a reversal of the signs of all the stresses is assumed to produce a reversal 
of sign of all the strain increments. This implies that C;; must be an odd function 
of the stress components and thus will vanish when all the s;; vanish. 

The material is supposed to become orthotropic under the stress o,; in the sense 
that the C,; can be represented as a power series in the stress deviator s;; with scalar 
coefficients. These coefficients are either constants or else functions of the invariants 
of s;;, i.e., functions of Jz and J3. It is convenient at this point to change from the sub- 
script notation for tensors to Gibbs’ notation; the tensor C,;; will be denoted by C 
and s;;by S. The multiplications indicated below are the usual matrix multiplications. 
Under the assumptions stated above, the tensor C can be written as 


C = >> aenss(J2, Js)S2**1. (25) 


n=0 


We note that only odd powers appear in Eq. (25) since C is assumed to be an odd 
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function of the stresses. Equation (25) can be simplified further by the Hamilton- 
Cayley theorem which states that the tensor S must satisfy its own characteristic 
equation.’ For the stress deviator S, this implies that 


S* = JS + Jal, (26) 


where I is the unit tensor. Through Eq. (26), we can reduce‘ any power of S greater 
than the second to a linear combination of I, S, and S? with coefficients which are 
functions of Jzand J3. For example, consider the reduction of the power S*. According 


to Eq. (26), 
S‘ = JS? + J; 


thus 
S'= 1,84 J8' =JS + IS + Jd. 
In general, we can rewrite Eq. (25) as 
C = a(J2, J3)S*? + b(J2, J3)S + c(Jo, J3)Jsl. 
We recall that C;;=0; since S is a deviator, this implies that 
2a(Jo, Js)J2 + 3c(J2, Js)J3 = 9, 


or 
(Js, Js) ile Ke 
C\J a, Sse eee Gg 2) ° 
‘ 3Js ; 


Consequently, 
C = a(Jo, Js) [S? — 3Je1] + b(J2, Js)S. 


The expression appearing in square brackets is just the tensor t;; which was defined 
in Eq. (13). Returning now to the subscript notation we can write the tensor C;; as 
Ci; = a(Jo, J 3)ti; + b(J2, J 3)Sij. 

A further simplification can be made by noting that C;; must be an odd function of the 

stress components. Since Jz is even, J; odd, t;; even, and s;; odd, we must have 
2. 2. 
Ci; ” p(J2, J3)Sij + q(Je, J 3)J stij- 
Thus the complete differential stress-strain relations can be written in the form$ 
de;;=(14+v)ds;;+(1 —2v)ds6;;+ [p(J2, Js)sij+q(J2, Js)Jatis]dJ2 when dJ22=0 ! an 
de; j= (1 +y)ds;;+ (1 — 2v)ds6; ;, when dJ2=0. 


5. Further study of the stress strain relations. Experimental determination. In 
this section we shall discuss the relation between the differential form of the stress- 


strain relation (cf. (27)) 





3M. Bocher, Introduction to higher algebra, The Macmillan Co., New York, 1907, p. 296. 

‘ This technique has been used recently by Marcus Reiner, Am. J. Math. 67, 350-362 (1945) and 
W. Prager, loc. cit. 

5 These relations contain, as special cases, the stress-strain laws developed by W. Prager, Proc. Fifth 
International Congress of Applied Mechanics, Cambridge, Mass., 1938, pp. 234-237, and by J. H. Laning 
in an unpublished paper (1942). 
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dei, = | p(Jn, Ia)8ij + QI, Js)Jatiz}dJ2, Jz 20 (28) 


and the integral form 
cj = F(x, Js) [PU Ja)sis + Qs Js)Jabis), (14) 
which holds only when the ratios of the principal stresses are kept constant during the 
loading process, i.e., if 
sig = RSs; (29) 
is fixed while k is the scalar variable. We shall then show how a series of 


where Sj; 
tests necessary to establish the Lode diagram will be sufficient to determine the stress- 


strain relations completely. First of all, it is convenient to bring out the homogeneity 
properties in the relations (28) and (14) by introducing the symbols 


2,3 2 
a =J3/Jo, Vij = J ti;/Jo, (30) 


where a@ is dimensionless, while y;; has the same dimensions as J2. The relation (14) 
can be written in the form 
Ad { ) . 
ij = A(J2, a) Siz + Bla)yiss, (31) 
where 


\(Jo, «) = F(J2, J3)P(Js, Ja), B(a) = J20(J2, Js)/P(J2, Js). (32) 


Note that 8 is independent of Jz, because of the homogeneity relation between P and 


Q established in Section 2. 
With a similar change of notation, the relation (28) can be written as 


de;; = G(Jo, a) {si; + B’(a)yis}dJo, dJ220 (33) 


where 
G(Jo, a) = p(J2, Js), B'(a) = Joq(Je, Js)/Pp(Ia, J). (34) 


Since we did not establish a homogeneity relation between p and q, we cannot immedi- 
lately conclude that 6’ is independent of J2. However, we shall see immediately that 


this is true and that indeed 
p’ =8. (35) 
We shall also show that G(Je, a) may be obtained from A(J2, a) by the relation 
r Or 
G(J2, a) = —+— : (36) 
22 OJ 


The relations (35) and (36) will then determine the differential relation (33) com- 
pletely once the integral relation (31) is known by a series of experiments of the 
special type (29). It is to be noted that the functions G(J2, aw) and B(a) in (33) deter- 
mined through the use of (31) will by no means restrict the application of (33) to 
processes connected in any manner with (29). 

To establish the relations (35) and (36), consider the application-of (31) and (33) 
to a process of the type (29). Let de/ be the change in ¢€// corresponding to a change dk. 
Then 
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dk 
dJ_ = 2): >? da = 0; (37) 
and (31) gives 
1 Or dk 
des; = — dJ2{ siz + Bris} + Al sis + Bris} —) 
aS 2 k 
while (33) yields 
dei; = GJ2, a) {sij + B’vis}dJ 2. 


Equating coefficients of s;; and y.;, we obtain the relations (35) and (36). 
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Fic. 2a. The a—8 diagram (Eq. (38)) of the experimental results of Taylor and Quinney for copper, 
aluminum, mild steel and decarburized mild steel. The data for mild steel are too scattered fora definite 


curve to be drawn. 


The experimental determination of the stress strain relations can then be reduced 
to that of (31) alone. This can be done by a series of tests of the type (29), which is 
of the class described by Lode,* Taylor-Quinney’ and Hohenemser-Prager.* Indeed, 








6 W. Lodge, Forschungsarbeiten a.d. Gebiete d. Ingenieurwesens, No. 303, VDI-Verlag, Berlin, 1928. 
7G. I. Taylor and H. Quinney, Phil. Trans. Roy. Soc. London (A) 230, 323-362 (1931). 

® K. Hohenemser and W. Prager, Z. angew. Math. Mech. 12, 1-14 (1932). An English translation of 
this paper is available as R.T.P. Translation No. 2468 (Durand Reprinting Committee, in care of Cali- 
fornia Institute of Technology, Pasadena 4, Calif.). 
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the relation 6(a) is merely another presentation of Lode’s diagram. It can be easily 
verified that a, 8 are related to Lode’s parameters® yp and »v by the relations 


4 p(9 — »*)* 
27 (3 + u?)* 
9(3 + y*)? 1 — v/p 

20—n) 4b] 23 








’ 


a= 


(38) 








B= 


This new system has the advantage that 6 gives directly the extent of deviation from 
“yon Mises’ second hypothesis” discussed by Taylor and Quinney, which is equiva- 
lent to putting 6=0. Indeed, one principal aim of Taylor and Quinney is to find out 
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Fic. 2b. The a—8 diagram (Eq. (38)) of the experimental results of Taylor and Quinney for lead, 
cadmium and glass. The data are too few to allow any curve to be drawn. 


this extent and is therefore to determine the value of 8. Figs. 2(a) and 2(b) show the 
results of Taylor and Quinney converted into the (a, 8) diagram. This diagram re- 
veals any experimental error more strongly, since 8 is essentially related to the slope 
of the (u, v) curve; e.g., 





9 W. Lode, loc. cit., pp. 1 and 12. 
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dv 9— 48 

— = —— for pl, 
du 9+ 28 

dy 


2p 
—=1+— for p=0. 
du 3 
It should be noted in passing that a must satisfy the inequality 


0 <a < 4/27 (39) 


to insure real values of p. 
Having determined B(a) from (38), we may determine A(J2, a) by noting that 
(cf. (31)) 


I, = heijei; = N22 {1 + 308 + 4a°B*}. (40) 


For each loading process given by (29) the value of a is fixed, and (40) gives the de- 
pendence of \? on Je if Iz is determined for given values of Je. A series of tests with 
different principal axes will then give the further dependence of \ on a. 

6. Concluding remarks. In closing, we note some of the limitations of the stress- 
strain relations developed in this paper. It has been pointed out previously that these 
equations have been developed to cover the case of one loading followed by at most 
one unloading. This restriction is quite essential, for relations (27) are not applicable 
for a second loading. For example, if we consider a simple tensile test, the stress-strain 
diagram obtained from (27) for the second loading would be a mere translation of the 
diagram for the first loading. This does not agree with the experimental results. Sec- 
ondly, we note that these equations apply only to metals which exhibit strain-harden- 
ing. They are not applicable, for example, to materials which yield under constant 
shearing stress'or satisfy von Mises’ yield condition, J2=const. 

It is hoped that the results presented here will provoke experimental work to 
test their validity. Among the various features which should be tested are two as- 
sumptions made in developing the differential stress-strain laws. The first hypothesis 
(cf. Section 4) states that the increments in strain are uniquely determined by the 
components of the stress tensor ¢;; and the increments do;; without reference to the 
previous history of loading provided only one loading has taken place followed by 
at most one unloading. The range in which such a hypothesis is valid should be ex- 
plored empirically. Secondly, the assumption involved in the transition from Eqs. 
(24) to (25) should be examined carefully. According to these two relations, the prin- 
cipal axes of the increment in permanent strain de{/ will coincide with the principal 
axes of the existing state of stress s;; independent of the increments in stress do;;, 
provided only loading takes place. This conclusion should be tested by experiment. 
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